Merge pull request #12 from spacetelescope/wcs_bug

WCS bugfix

CHANGES.rst

@@ -2,6 +2,7 @@
 ----------------
 
 - Adding more unit tests and coveralls setup [#11]
+- Adding workaround for FFIs with bad WCS info [#12]
 
 
 0.3 (2019-05-03)

astrocut/cube_cut.py

@@ -72,7 +72,7 @@ def __init__(self):
                          "VIGNAPP": [None, "vignetting or collimator correction applied"]}
 
 
-    def _parse_table_info(self, table_header, table_row):
+    def _parse_table_info(self, table_header, table_data, verbose=False):
         """
         Takes the header and one entry from the cube table of image header data,
         builds a WCS object that encalpsulates the given WCS information,
@@ -88,7 +88,23 @@ def _parse_table_info(self, table_header, table_row):
             One row from the cube image header data table.
         """
 
-        # Turning the table row inot a new header object
+        data_ind = int(len(table_data)/2)  # using the middle file for table info
+        table_row = None
+
+        # Making sure we have a row with wcs info
+        while table_row is None:
+            table_row = table_data[data_ind]
+            ra_col = int([x for x in table_header.cards if x[1] == "CTYPE1"][0][0].replace("TTYPE", "")) - 1
+            if "RA" not in table_row[ra_col]:
+                table_row is None
+                data_ind += 1
+                if data_ind == len(table_data):
+                    raise wcs.NoWcsKeywordsFoundError("No FFI rows contain valid WCS keywords.")
+
+        if verbose:
+            print("Using WCS from row {} out of {}".format(data_ind, len(table_data)))
+
+        # Turning the table row into a new header object
         wcs_header = fits.header.Header()
 
         for header_key, header_val in table_header.items():
@@ -97,6 +113,11 @@ def _parse_table_info(self, table_header, table_row):
 
             col_num = int(header_key[5:]) - 1
             tform = table_header[header_key.replace("TTYPE", "TFORM")]
+
+            wcs_val = table_row[col_num]
+            if (not isinstance(wcs_val, str)) and (np.isnan(wcs_val)):
+                continue  # Just skip nans
+
             if "A" in tform:
                 wcs_header[header_val] = str(table_row[col_num])
             elif "D" in tform:
@@ -135,7 +156,10 @@ def _get_cutout_limits(self, cutout_size):
         """
 
         # Note: This is returning the center pixel in 1-up
-        center_pixel = self.center_coord.to_pixel(self.cube_wcs, 1)
+        try:
+            center_pixel = self.center_coord.to_pixel(self.cube_wcs, 1)
+        except wcs.NoConvergence:  # If wcs can't converge, center coordinate is far from the footprint
+            raise InvalidQueryError("Cutout location is not in cube footprint!")
 
         lims = np.zeros((2, 2), dtype=int)
 
@@ -619,8 +643,7 @@ def cube_cut(self, cube_file, coordinates, cutout_size,
         cube = fits.open(cube_file) 
 
         # Get the info we need from the data table
-        data_ind = int(len(cube[2].data)/2)  # using the middle file for table info
-        self._parse_table_info(cube[2].header, cube[2].data[data_ind])
+        self._parse_table_info(cube[2].header, cube[2].data, verbose)
 
         if isinstance(coordinates, SkyCoord):
             self.center_coord = coordinates

astrocut/make_cube.py

@@ -93,56 +93,79 @@ def make_cube(self, file_list, cube_file="img-cube.fits", sector=None, verbose=T
             startTime = time()
 
         # Getting the time sorted indices for the files
+        # and a good header to build the table colums (needs to have WCS keywords)
         file_list = np.array(file_list)
         start_times = np.zeros(len(file_list))
+        good_header_ind = None
 
         for i, ffi in enumerate(file_list):
-            start_times[i] = fits.getval(ffi, 'TSTART', 1)  # TODO: optionally pass this in?
+            ffi_data = fits.open(ffi)
+
+            start_times[i] = ffi_data[1].header.get("TSTART")  # TODO: optionally pass this in?
+
+            if good_header_ind is None:  # Only check this if we don't already have it
+                if ffi_data[1].header.get("CTYPE1"):  # Checking for WCS info
+                    good_header_ind = i
+
+            ffi_data.close()
+
+
 
         sorted_indices = np.argsort(start_times)      
 
-        # these will be the arrays and headers
+        # These will be the arrays and headers
         img_cube = None
         img_info_table = None
         primary_header = None
-        secondary_header = None
+
+        # Setting up the image info table
+        fle = file_list[good_header_ind]
+
+        if verbose:
+            print("Using {} to initialize the image header table.".format(path.basename(fle)))
+
+        ffi_data = fits.open(fle)
+
+        # The image specific header information will be saved in a table in the second extension
+        secondary_header = ffi_data[1].header
+
+        # set up the image info table
+        cols = []
+        for kwd, val, cmt in secondary_header.cards: 
+            if type(val) == str:
+                tpe = "S" + str(len(val))  # TODO: Maybe switch to U?
+            elif type(val) == int:
+                tpe = np.int32
+            else:
+                tpe = np.float32
+
+            cols.append(Column(name=kwd, dtype=tpe, length=len(file_list), meta={"comment": cmt}))
+
+        cols.append(Column(name="FFI_FILE", dtype="S" + str(len(path.basename(fle))), length=len(file_list)))
+
+        img_info_table = Table(cols)
+
+        ffi_data.close()
+
 
         # Loop through files
         for i, fle in enumerate(file_list[sorted_indices]):
 
             ffi_data = fits.open(fle)
 
-            # if the arrays/headers aren't initialized do it now
+            # if the arrays/header aren't initialized do it now
             if img_cube is None:
 
                 # We use the primary header from the first file as the cube primary header
                 # and will add in information about the time of the final observation at the end
                 primary_header = self._make_primary_header(ffi_data[0].header, ffi_data[1].header, sector=sector)
 
-                # The image specific header information will be saved in a table in the second extension
-                secondary_header = ffi_data[1].header
-
                 ffi_img = ffi_data[1].data
 
                 # set up the cube array
                 img_cube = np.full((ffi_img.shape[0], ffi_img.shape[1], len(file_list), 2),
                                    np.nan, dtype=np.float32)
 
-                # set up the image info table
-                cols = []
-                for kwd, val, cmt in secondary_header.cards: 
-                    if type(val) == str:
-                        tpe = "S" + str(len(val))  # TODO: Maybe switch to U?
-                    elif type(val) == int:
-                        tpe = np.int32
-                    else:
-                        tpe = np.float32
-
-                    cols.append(Column(name=kwd, dtype=tpe, length=len(file_list), meta={"comment": cmt}))
-
-                cols.append(Column(name="FFI_FILE", dtype="S" + str(len(path.basename(fle))), length=len(file_list)))
-
-                img_info_table = Table(cols)
 
             # add the image and info to the arrays
             img_cube[:, :, i, 0] = ffi_data[1].data