[Moving Objects Telescope Exposure Crossmatch] Find which known minor-planet bodies lie within any telescope exposure with a square FOV #pyephem #moving_object #fov #atlas #crossmatch

minor-planet-square-exposure-crossmatch.py

#!/usr/local/bin/python
# encoding: utf-8
"""
*Find which known minor-planet bodies lie within any telescope exposure with a square FOV*

:Author:
    David Young

:Date Created:
    June 30, 2017

Usage:
    minor-planet-square-exposure-crossmatch <ra> <dec> <mjd> <width> <obsCode>

Options:
    ra                    right ascension of the exposure field centre
    dec                   declination of the exposure field centre
    mjd                   the mjd of the expsure (preferably the mid-exposure mjd)
    width                 the on-sky width of one side of the square exposure
    obsCode               the observatory code

    -h, --help            show this help message
    -v, --version         show version
    -s, --settings        the settings file
"""
################# GLOBAL IMPORTS ####################
import sys
import os
from fundamentals import tools
import codecs
import ephem
import math
from astrocalc.coords import unit_conversion, separations
import healpy as hp
import numpy as np


def main(arguments=None):
    """
    *The main function used when ``minor-planet-square-exposure-crossmatch.py`` is run as a single script from the cl*
    """

    # SETUP THE COMMAND-LINE UTIL SETTINGS
    su = tools(
        arguments=arguments,
        docString=__doc__,
        logLevel="WARNING",
        options_first=True,
        projectName=False
    )
    arguments, settings, log, dbConn = su.setup()

    # UNPACK REMAINING CL ARGUMENTS USING `EXEC` TO SETUP THE VARIABLE NAMES
    # AUTOMATICALLY
    for arg, val in arguments.iteritems():
        if arg[0] == "-":
            varname = arg.replace("-", "") + "Flag"
        else:
            varname = arg.replace("<", "").replace(">", "")
        if isinstance(val, str) or isinstance(val, unicode):
            exec(varname + " = '%s'" % (val,))
        else:
            exec(varname + " = %s" % (val,))
        if arg == "--dbConn":
            dbConn = val
        log.debug('%s = %s' % (varname, val,))

    # MAKE SURE HEALPIX SMALL ENOUGH TO MATCH FOOTPRINTS CORRECTLY
    nside = 1024
    pi = (4 * math.atan(1.0))
    DEG_TO_RAD_FACTOR = pi / 180.0
    RAD_TO_DEG_FACTOR = 180.0 / pi

    obs = ephem.Observer()
    obs.lat = obsLat
    obs.lon = obsLon
    # obs.elevation = 3000
    # obs.epoch = ephem.J2000
    mjd = float(mjd)
    raFc = float(ra)
    decFc = float(dec)
    tileSide = float(width)

    # READ IN MPC EPHEMERIS DATABASE
    pathToReadFile = "MPCORB.edb"
    try:
        log.debug("attempting to open the file %s" % (pathToReadFile,))
        readFile = codecs.open(pathToReadFile, encoding='utf-8', mode='r')
        edb = readFile.readlines()
        readFile.close()
    except IOError, e:
        message = 'could not open the file %s' % (pathToReadFile,)
        log.critical(message)
        raise IOError(message)
    readFile.close()

    lastDes = 0
    numbered = True
    ra = []
    dec = []
    healpix = []
    objects = []

    # GENERATE THE POSITIONS FOR EACH SOURCE IN EPH DB FILE
    for line in edb:
        if line[0] == "#":
            continue

        # PyEphem works in DUBLIN JD, to convert from MJD subtract 15019.5
        minorPlanet = ephem.readdb(line.strip())
        obs.date = mjd - 15019.5
        minorPlanet.compute(obs)

        nameList = minorPlanet.name.split()
        if numbered == True and int(nameList[0]) > lastDes:
            lastDes = int(nameList[0])
            desig = int(nameList[0])
        else:
            numbered = False
            desig = minorPlanet.name

        # THE FOLLOW ARE SET FOR ALL BODIES
        # Astrometric geocentric right ascension for the epoch specified
        objects.append(minorPlanet.name)
        if "276409" in str(desig):
            print minorPlanet.ra, minorPlanet.dec
            print float(minorPlanet.ra), float(minorPlanet.dec)
            print obs.date
            print minorPlanet.ra * RAD_TO_DEG_FACTOR, minorPlanet.dec * RAD_TO_DEG_FACTOR
            print minorPlanet.a_ra * RAD_TO_DEG_FACTOR, minorPlanet.a_dec * RAD_TO_DEG_FACTOR
        ra.append(minorPlanet.ra * RAD_TO_DEG_FACTOR)
        dec.append(minorPlanet.dec * RAD_TO_DEG_FACTOR)

        # Astrometric geocentric declination for the epoch specified
        # print minorPlanet.a_dec * RAD_TO_DEG_FACTOR
        # print minorPlanet.ra  # Apparent geocentric right ascension for the epoch-of-date
        # print minorPlanet.g_dec  # Apparent geocentric declination for the epoch-of-date
        # print minorPlanet.elong  # Elongation (angle to sun)
        # print minorPlanet.mag  # Magnitude
        # print minorPlanet.size  # Size (diameter in arcseconds)
        # print minorPlanet.radius  # Size (radius as an angle)
        # # print minorPlanet.circuminorPlanetolar  # whether it stays above the horizon  -- requires observer not date
        # # print minorPlanet.neverup  # whether is stays below the horizon  -- requires
        # # observer not date

        # print minorPlanet.hlon  # Heliocentric longitude (see next paragraph)
        # print minorPlanet.hlat  # Heliocentric latitude (see next paragraph)
        # print minorPlanet.sun_distance  # Distance to Sun (AU)
        # print minorPlanet.earth_distance  # Distance to Earth (AU)
        # print minorPlanet.phase  # Percent of surface illuminated

    # GERNEATE THE SKY-LOCATIONS AS A HEALPIXEL ID
    ra = np.array(ra)
    dec = np.array(dec)
    healpix = hp.ang2pix(nside, theta=ra, phi=dec, lonlat=True)

    # GENERATE THE EXPOSURE HEALPIX ID MAP
    decCorners = (decFc - tileSide / 2,
                  decFc + tileSide / 2)
    corners = []
    for d in decCorners:
        if d > 90.:
            d = 180. - d
        elif d < -90.:
            d = -180 - d
        raCorners = (raFc - (tileSide / 2) / np.cos(d * DEG_TO_RAD_FACTOR),
                     raFc + (tileSide / 2) / np.cos(d * DEG_TO_RAD_FACTOR))
        for r in raCorners:
            if r > 360.:
                r = 720. - r
            elif r < 0.:
                r = 360. + r
            corners.append(hp.ang2vec(r, d, lonlat=True))

    # NEAR THE POLES RETURN SQUARE INTO TRIANGE - ALMOST DEGENERATE
    pole = False
    for d in decCorners:
        if d > 87.0 or d < -87.0:
            pole = True

    if pole == True:
        corners = corners[1:]
    else:
        # FLIP CORNERS 3 & 4 SO HEALPY UNDERSTANDS POLYGON SHAPE
        corners = [corners[0], corners[1],
                   corners[3], corners[2]]

    # RETURN HEALPIXELS IN EXPOSURE AREA
    expPixels = hp.query_polygon(nside, np.array(
        corners))

    from collections import defaultdict
    dicto = defaultdict(list)  # dictionary stores all the relevant coordinates
    # so you don't have to search for them later

    for ind, (p, r, d, o) in enumerate(zip(healpix, ra, dec, objects)):
        dicto[p].append([o, r, d])

    matches = []
    for ele in expPixels:
        if ele in dicto:
            matches.append(dicto[ele])

    pathToWriteFile = "dryx.kastx"
    try:
        log.debug("attempting to open the file %s" % (pathToWriteFile,))
        writeFile = codecs.open(pathToWriteFile, encoding='utf-8', mode='w')
    except IOError, e:
        message = 'could not open the file %s' % (pathToWriteFile,)
        log.critical(message)
        raise IOError(message)

    writeFile.write("# OBJECT, RA, DEC\n")

    for m in matches:
        o = m[0][0]
        r = m[0][1]
        d = m[0][2]
        writeFile.write("""%(o)s, %(r)0.5f, %(d)0.5f\n""" % locals())

    writeFile.close()

    return


if __name__ == '__main__':
    main()