Merge branch 'skyfit2fix' into dev

RMS/Astrometry/ApplyAstrometry.py

@@ -41,7 +41,7 @@
 from RMS.Formats.FTPdetectinfo import readFTPdetectinfo, writeFTPdetectinfo
 from RMS.Formats.FFfile import filenameToDatetime
 import RMS.Formats.Platepar
-from RMS.Math import angularSeparation
+from RMS.Math import angularSeparation, polarToCartesian, cartesianToPolar
 
 # Import Cython functions
 import pyximport
@@ -293,7 +293,6 @@ def photometryFitRobust(px_intens_list, radius_list, catalog_mags, fixed_vignett
 
 
 
-
 def computeFOVSize(platepar):
     """ Computes the size of the FOV in deg from the given platepar.
 
@@ -305,30 +304,70 @@ def computeFOVSize(platepar):
     """
 
     # Construct poinits on the middle of every side of the image
-    time_data = np.array(4*[jd2Date(platepar.JD)])
-    x_data = np.array([0, platepar.X_res, platepar.X_res/2, platepar.X_res/2])
-    y_data = np.array([platepar.Y_res/2, platepar.Y_res/2, 0, platepar.Y_res])
-    level_data = np.ones(4)
+    x_data = np.array([               0,  platepar.X_res,  platepar.X_res/2, platepar.X_res/2, platepar.X_res/2.0])
+    y_data = np.array([platepar.Y_res/2, platepar.Y_res/2,                0, platepar.Y_res,   platepar.Y_res/2.0])
+    time_data = np.array(len(x_data)*[jd2Date(platepar.JD)])
+    level_data = np.ones(len(x_data))
 
     # Compute RA/Dec of the points
-    _, ra_data, dec_data, _ = xyToRaDecPP(time_data, x_data, y_data, level_data, platepar)
+    _, ra_data, dec_data, _ = xyToRaDecPP(time_data, x_data, y_data, level_data, platepar, \
+        extinction_correction=False)
 
-    ra1, ra2, ra3, ra4 = ra_data
-    dec1, dec2, dec3, dec4 = dec_data
+    ra1, ra2, ra3, ra4, ra_mid = ra_data
+    dec1, dec2, dec3, dec4, dec_mid = dec_data
 
     # Compute horizontal FOV
-    fov_h = np.degrees(angularSeparation(np.radians(ra1), np.radians(dec1), np.radians(ra2), \
-        np.radians(dec2)))
+    fov_hl = np.degrees(angularSeparation(np.radians(ra1), np.radians(dec1), np.radians(ra_mid), \
+        np.radians(dec_mid)))
+    fov_hr = np.degrees(angularSeparation(np.radians(ra2), np.radians(dec2), np.radians(ra_mid), \
+        np.radians(dec_mid)))
+    fov_h = fov_hl + fov_hr
 
     # Compute vertical FOV
-    fov_v = np.degrees(angularSeparation(np.radians(ra3), np.radians(dec3), np.radians(ra4), \
-        np.radians(dec4)))
-
+    fov_vu = np.degrees(angularSeparation(np.radians(ra3), np.radians(dec3), np.radians(ra_mid), \
+        np.radians(dec_mid)))
+    fov_vd = np.degrees(angularSeparation(np.radians(ra4), np.radians(dec4), np.radians(ra_mid), \
+        np.radians(dec_mid)))
+    fov_v = fov_vu + fov_vd
 
     return fov_h, fov_v
 
 
 
+def getFOVSelectionRadius(platepar):
+    """ Get a radius around the centre of the FOV which includes the FOV, but excludes stars outside the FOV.
+    Arguments:
+        platepar: [Platepar instance]
+
+    Return:
+        fov_radius: [float] Radius in degrees.
+    """
+
+    # Construct poinits on the middle of every side of the image
+    x_data = np.array([0, platepar.X_res, platepar.X_res,              0, platepar.X_res/2.0])
+    y_data = np.array([0, platepar.Y_res,              0, platepar.Y_res, platepar.Y_res/2.0])
+    time_data = np.array(len(x_data)*[jd2Date(platepar.JD)])
+    level_data = np.ones(len(x_data))
+
+    # Compute RA/Dec of the points
+    _, ra_data, dec_data, _ = xyToRaDecPP(time_data, x_data, y_data, level_data, platepar, \
+        extinction_correction=False)
+
+    ra1, ra2, ra3, ra4, ra_mid = ra_data
+    dec1, dec2, dec3, dec4, dec_mid = dec_data
+
+    # Angular separation between the centre of the FOV and corners
+    ul_sep = np.degrees(angularSeparation(np.radians(ra1), np.radians(dec1), np.radians(ra_mid), np.radians(dec_mid)))
+    lr_sep = np.degrees(angularSeparation(np.radians(ra2), np.radians(dec2), np.radians(ra_mid), np.radians(dec_mid)))
+    ur_sep = np.degrees(angularSeparation(np.radians(ra3), np.radians(dec3), np.radians(ra_mid), np.radians(dec_mid)))
+    ll_sep = np.degrees(angularSeparation(np.radians(ra4), np.radians(dec4), np.radians(ra_mid), np.radians(dec_mid)))
+
+    # Take the average radius
+    fov_radius = np.mean([ul_sep, lr_sep, ur_sep, ll_sep])
+
+    return fov_radius
+
+
 
 def rotationWrtHorizon(platepar):
     """ Given the platepar, compute the rotation of the FOV with respect to the horizon.

RMS/Astrometry/AstrometryNetNova.py

@@ -9,6 +9,7 @@
 
 import os
 import sys
+import copy
 import time
 import base64
 import json
@@ -338,6 +339,16 @@ def novaAstrometryNetSolve(ff_file_path=None, img=None, x_data=None, y_data=None
     """
 
 
+    def _printWebLink(stat, first_status=None):
+
+        if first_status is not None:
+            if not len(stat.get("user_images", "")):
+                stat = first_status
+
+        if len(stat.get("user_images", "")):
+            print("Link to web page: http://nova.astrometry.net/user_images/{:d}".format(stat.get("user_images", "")[0]))
+
+
     # Read the FF file, if given
     if ff_file_path is not None:
 
@@ -355,7 +366,7 @@ def novaAstrometryNetSolve(ff_file_path=None, img=None, x_data=None, y_data=None
 
         # Save the avepixel as a memory file
         file_handle = BytesIO()
-        pil_img = Image.fromarray(img)
+        pil_img = Image.fromarray(img.T)
 
         # Save image to memory as JPG
         pil_img.save(file_handle, format='JPEG')
@@ -376,7 +387,7 @@ def novaAstrometryNetSolve(ff_file_path=None, img=None, x_data=None, y_data=None
 
     # Add keyword arguments
     kwargs = {}
-    kwargs['publicly_visible'] = 'n'
+    kwargs['publicly_visible'] = 'y'
     kwargs['crpix_center'] = True
     kwargs['tweak_order'] = 3
 
@@ -418,8 +429,9 @@ def novaAstrometryNetSolve(ff_file_path=None, img=None, x_data=None, y_data=None
     tries = 0
     while True:
 
-        # Limit the number of checking if the fiels is solved, so the script does not get stuck
+        # Limit the number of checking if the field is solved, so the script does not get stuck
         if tries > solution_tries:
+            _printWebLink(stat)
             return None
 
         stat = c.sub_status(sub_id, justdict=True)
@@ -441,6 +453,9 @@ def novaAstrometryNetSolve(ff_file_path=None, img=None, x_data=None, y_data=None
 
         tries += 1
 
+
+    first_status = copy.deepcopy(stat)
+
     # Get results
     get_results_tries = 10
     get_solution_tries = 30
@@ -451,10 +466,12 @@ def novaAstrometryNetSolve(ff_file_path=None, img=None, x_data=None, y_data=None
         # Limit the number of tries of getting the results, so the script does not get stuck
         if results_tries > get_results_tries:
             print('Too many tries in getting the results!')
+            _printWebLink(stat, first_status=first_status)
             return None
 
         if solution_tries > get_solution_tries:
             print('Waiting too long for the solution!')
+            _printWebLink(stat, first_status=first_status)
             return None
 
         # Get the job status
@@ -470,6 +487,9 @@ def novaAstrometryNetSolve(ff_file_path=None, img=None, x_data=None, y_data=None
 
         elif stat.get('status','') in ['failure']:
             print('Failed to find a solution!')
+
+            _printWebLink(stat, first_status=first_status)
+
             return None
 
         # Wait until the job is solved

RMS/Astrometry/CheckFit.py

@@ -20,7 +20,7 @@
 from RMS.Formats import CALSTARS
 from RMS.Formats import StarCatalog
 from RMS.Formats import FFfile
-from RMS.Astrometry.ApplyAstrometry import raDecToXYPP, xyToRaDecPP, rotationWrtHorizon
+from RMS.Astrometry.ApplyAstrometry import raDecToXYPP, xyToRaDecPP, rotationWrtHorizon, getFOVSelectionRadius
 from RMS.Astrometry.Conversions import date2JD, jd2Date, raDec2AltAz
 from RMS.Astrometry.FFTalign import alignPlatepar
 from RMS.Math import angularSeparation
@@ -78,15 +78,7 @@ def matchStarsResiduals(config, platepar, catalog_stars, star_dict, match_radius
 
 
     # Estimate the FOV radius
-    fov_w = platepar.X_res/platepar.F_scale
-    fov_h = platepar.Y_res/platepar.F_scale
-
-    fov_radius = np.sqrt((fov_w/2)**2 + (fov_h/2)**2)
-
-    # print('fscale', platepar.F_scale)
-    # print('FOV w:', fov_w)
-    # print('FOV h:', fov_h)
-    # print('FOV radius:', fov_radius)
+    fov_radius = getFOVSelectionRadius(platepar)
 
 
     # Dictionary containing the matched stars, the keys are JDs of every image

RMS/Astrometry/FFTalign.py

@@ -238,7 +238,7 @@ def alignPlatepar(config, platepar, calstars_time, calstars_coords, scale_update
 
         # Calculate the FOV radius in degrees
         fov_y, fov_x = ApplyAstrometry.computeFOVSize(platepar)
-        fov_radius = np.sqrt(fov_x**2 + fov_y**2)
+        fov_radius = ApplyAstrometry.getFOVSelectionRadius(platepar)
 
         # Take only those stars which are inside the FOV
         filtered_indices, _ = subsetCatalog(catalog_stars, ra_centre, dec_centre, jd, platepar.lat, \

RMS/Astrometry/SkyFit.py

@@ -52,7 +52,8 @@
 
 from RMS.Astrometry.ApplyAstrometry import xyToRaDecPP, raDecToXYPP, \
     rotationWrtHorizon, rotationWrtHorizonToPosAngle, computeFOVSize, photomLine, photometryFit, \
-    rotationWrtStandard, rotationWrtStandardToPosAngle, correctVignetting, extinctionCorrectionTrueToApparent
+    rotationWrtStandard, rotationWrtStandardToPosAngle, correctVignetting, extinctionCorrectionTrueToApparent,
+    getFOVSelectionRadius
 from RMS.Astrometry.AstrometryNetNova import novaAstrometryNetSolve
 from RMS.Astrometry.Conversions import date2JD, JD2HourAngle
 import RMS.ConfigReader as cr
@@ -2222,8 +2223,7 @@ def filterCatalogStarsInsideFOV(self, catalog_stars):
         ra_centre, dec_centre = self.computeCentreRADec()
 
         # Calculate the FOV radius in degrees
-        fov_y, fov_x = computeFOVSize(self.platepar)
-        fov_radius = np.sqrt(fov_x**2 + fov_y**2)
+        fov_radius = getFOVSelectionRadius(self.platepar)
 
         # Compute the current Julian date
         jd = date2JD(*self.img_handle.currentTime())