Merge pull request #6 from stevenstetzler/arcsat_fixes

Support ARCSAT observatory images and add console script

setup.py

@@ -1,3 +1,9 @@
 from setuptools import setup
 
-setup()
+setup(
+    entry_points={
+        "console_scripts": [
+            "koffi = koffi.script:main"
+        ]
+    }
+)

src/koffi/koffi.py

@@ -30,7 +30,6 @@
 from koffi_tools.image_metadata import *
 from koffi_tools.potential_source import *
 
-
 def skybot_search_frame(image):
     """
     Gets all known objects within the frame of a single FITS image from SkyBoT.
@@ -162,7 +161,7 @@ def create_jpl_query_string(image):
 
     # Create a string of data for the observatory.
     if image.obs_code:
-        obs_str = "mpc-code=%s" % self.obs_code
+        obs_str = "mpc-code=%s" % image.obs_code
     else:
         obs_str = "lat=%f&lon=%f&alt=%f" % (image.obs_lat, image.obs_long, image.obs_alt)
 
@@ -189,7 +188,7 @@ def create_jpl_query_string(image):
         dec_str = "fov-dec-lim=M%02i-%02i-%05.2f" % (-dec_dms_L[0], -dec_dms_L[1], -dec_dms_L[2])
     dec_dms_H = Angle(image.center.dec + image.dec_radius()).dms
     if dec_dms_H[0] >= 0:
-        dec_str = "%s,02i-%02i-%05.2f" % (dec_str, dec_dms_H[0], dec_dms_H[1], dec_dms_H[2])
+        dec_str = "%s,%02i-%02i-%05.2f" % (dec_str, dec_dms_H[0], dec_dms_H[1], dec_dms_H[2])
     else:
         dec_str = "%s,M%02i-%02i-%05.2f" % (dec_str, -dec_dms_H[0], -dec_dms_H[1], -dec_dms_H[2])
 
@@ -223,6 +222,9 @@ def jpl_search_frame(image):
         feed = url.read().decode("utf-8")
         results = json.loads(feed)
 
+        if 'warning' in results.keys() and results['warning'] == "no matching records":
+            return []
+
         num_results = results["n_second_pass"]
         for item in results["data_second_pass"]:
             name = item[0]

src/koffi/script.py

@@ -0,0 +1,48 @@
+from .koffi import ImageMetadata, skybot_search_frame, jpl_search_frame
+from astropy.table import QTable
+from astropy import units as u
+import argparse
+import sys
+
+
+def main():
+    parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    parser.add_argument("filename", type=str, help="The fits filts to open")
+    parser.add_argument("--query", type=str, default="JPL", help="The service to query, JPL or SkyBot")
+    parser.add_argument("--format", type=str, default="QTable", help="The output format to use")
+
+    args = parser.parse_args()
+
+    image = ImageMetadata(args.filename)
+
+    if args.query == "JPL":
+        search = jpl_search_frame
+    elif args.query == "SkyBot":
+        search = skybot_search_frame
+
+    objects = search(image)
+    names = []
+    ras = []
+    decs = []
+    xs = []
+    ys = []
+    for name, coord in objects:
+        names.append(name)
+        ras.append(coord.ra.degree * u.degree)
+        decs.append(coord.dec.degree * u.degree)
+        x, y = image.wcs.world_to_pixel(coord)
+        xs.append(x)
+        ys.append(y)
+
+    table = QTable(
+        data=[names, ras, decs, xs, ys], 
+        names=["Object Name", "RA", "Dec", "x", "y"],
+    )
+    if args.format == "QTable":
+        print(table)
+    else:
+        table.write(sys.stdout, format=args.format)
+
+
+if __name__ == "__main__":
+    main()

src/koffi_tools/image_metadata.py

@@ -1,6 +1,7 @@
 from astropy.io import fits
-from astropy.time import Time
+from astropy.time import Time, TimeDelta
 from astropy.wcs import WCS
+from astroquery.mpc import MPC
 
 
 class ImageMetadata:
@@ -57,16 +58,19 @@ def populate_from_fits_file(self, filename, mjd_key="MJD_OBS", mjd_val=None):
                 # Patch for DECam data
                 if "DATE-AVG" in hdu_list[0].header:
                     self.set_epoch(Time(hdu_list[0].header["DATE-AVG"], format="isot"))
+                elif "DATE-OBS" in hdu_list[0].header and "EXPTIME" in hdu_list[0].header:
+                    self.set_epoch(Time(hdu_list[0].header["DATE-OBS"], format="isot") + TimeDelta(hdu_list[0].header["EXPTIME"], format="sec")/2)
                 elif self.get_header_element(mjd_key) is not None:
                     self.set_epoch(Time(self.get_header_element(mjd_key), format="mjd"))
 
             # Extract information about the location of the observatory.
             # Since this doesn't seem to be standardized, we try some
             # documented versions.
             observat = self.get_header_element("OBSERVAT")
+            observatories = MPC.get_observatory_codes() # get list of MPC-defined observatory codes
             obs_lat = self.get_header_element("OBS-LAT")
-            lat_obs = self.get_header_element("LAT_OBS")
-            if observat is not None:
+            lat_obs = self.get_header_element("LAT-OBS")
+            if observat is not None and observat in observatories['Code']:
                 self.obs_code = observat
                 self.obs_loc_set = True
             elif obs_lat is not None:
@@ -76,8 +80,8 @@ def populate_from_fits_file(self, filename, mjd_key="MJD_OBS", mjd_val=None):
                 self.obs_loc_set = True
             elif lat_obs is not None:
                 self.obs_lat = float(lat_obs)
-                self.obs_long = float(self.get_header_element("LONG_OBS"))
-                self.obs_alt = float(self.get_header_element("ALT_OBS"))
+                self.obs_long = float(self.get_header_element("LONG-OBS"))
+                self.obs_alt = float(self.get_header_element("ALT-OBS"))
                 self.obs_loc_set = True
             else:
                 self.obs_loc_set = False