Adding error test and fixes to MOSFIRE methods

potpyri/instruments/MOSFIRE.py

@@ -92,7 +92,7 @@ def get_saturation(self, hdr):
         return(hdr['SATURATE']*hdr['SYSGAIN'])
 
     def raw_format(self, proc):
-        return('MF.*.fits.gz')
+        return('*.fits.gz')
 
     def get_filter(self, hdr):
         filt = hdr['FILTER'].replace(' ','').split('_')[0]
@@ -117,7 +117,13 @@ def get_number(self, header):
 
     def import_image(self, filename, amp, log=None):
 
-        raw = CCDData.read(filename, unit=u.adu)
+        with fits.open(filename) as hdr:
+            header = hdr['SCI'].header
+            data = hdr['SCI'].data
+
+        del header['BUNIT']
+
+        raw = CCDData(data, header=header, unit=u.adu)
         red = ccdproc.ccd_process(raw, 
             gain=self.get_gain(raw.header)*u.electron/u.adu, 
             readnoise=self.get_rdnoise(raw.header)*u.electron)

potpyri/primitives/image_procs.py

@@ -747,8 +747,8 @@ def create_error(science_data, mask_data, rdnoise):
     mask = np.isnan(error)
     error[mask] = np.nanmedian(error)
     mask = error < 0.0
-    error = np.nanmedian(error)
-    maxval = np.max(img_data)
+    error[mask] = np.nanmedian(error)
+    maxval = np.float32(hdu[0].header['SATURATE'])
     mask = np.isinf(error)
     error[mask] = maxval
 

potpyri/primitives/solve_wcs.py

@@ -134,7 +134,12 @@ def solve_astrometry(file, tel, binn, paths, radius=0.5, replace=True,
 
     exten = '.'+file.split('.')[-1]
 
-    check_pairs = [('CRVAL1','CRVAL2'),('RA','DEC'),('OBJCTRA','OBJCTDEC')]
+    if tel.name.upper=='BINOSPEC':
+        check_pairs = [('CRVAL1','CRVAL2'),('RA','DEC'),('OBJCTRA','OBJCTDEC')]
+    else:
+        check_pairs = [('RA','DEC'),('OBJCTRA','OBJCTDEC'),('TARGRA','TARGDEC'),
+            ('CRVAL1','CRVAL2')]
+
     coord = None
 
     for pair in check_pairs:
@@ -198,7 +203,7 @@ def solve_astrometry(file, tel, binn, paths, radius=0.5, replace=True,
 
     tries = 1
     good = False
-    while tries < 5 and not good:
+    while tries < 6 and not good:
         input_args = f'{args} {extra_opts}'
 
         if log: 
@@ -230,6 +235,15 @@ def solve_astrometry(file, tel, binn, paths, radius=0.5, replace=True,
             elif tries==3:
                 extra_opts='--no-verify'
             elif tries==4:
+                # Try without source extractor
+                extra_opts='--no-verify'
+                args = '--scale-units arcsecperpix '
+                args += f'--scale-low {scale_low} --scale-high {scale_high} '
+                args += f'--no-plots -T '
+                args += f'--overwrite -N {newfile} --dir {directory} '
+                args += f'--ra {ra} --dec {dec} '
+                args += f' --radius {radius} '
+            elif tries==5:
                 # Try with no constraint on RA/Dec
                 args = '--scale-units arcsecperpix '
                 args += f'--scale-low {scale_low} --scale-high {scale_high} '

potpyri/primitives/sort_files.py

@@ -248,7 +248,7 @@ def sort_files(files, file_list, tel, paths, incl_bad=False, log=None):
             moved_path = paths['bad']
             if os.path.dirname(f)!=moved_path: shutil.move(f, paths['bad'])
             continue
-        except (TypeError, gzip.BadGzipFile, zlib.error):
+        except (TypeError, gzip.BadGzipFile, zlib.error, OSError):
             if log: 
                 log.error(f'Moving file {f} to bad due to corrupted data.')
             else:

tests/test_error.py

@@ -0,0 +1,54 @@
+from potpyri.utils import options
+from potpyri.utils import logger
+from potpyri.primitives import image_procs
+from potpyri.instruments import instrument_getter
+
+import numpy as np
+import os
+
+from ccdproc import CCDData
+
+from astropy.io import fits
+from astropy.wcs import WCS
+from astropy import units as u
+
+from . import utils
+
+def test_error(tmp_path):
+
+    instrument = 'MOSFIRE'
+    file_list_name = 'files.txt'
+
+    # Science file (just tellurics)
+    tel = instrument_getter(instrument)
+    paths = options.add_paths(tmp_path, file_list_name, tel)
+
+    # Generate log file in corresponding directory for log
+    log = logger.get_log(paths['log'])
+
+    aligned_data = []
+    masks = []
+    errors = []
+
+    rdnoise = 10.
+    imgval = 1000.
+
+    imdata = np.ones(shape=(2048, 2048), dtype=np.float32) * imgval
+    maskim = np.zeros(shape=(2048, 2048), dtype=np.uint16)
+    maskhdu = fits.PrimaryHDU(maskim)
+
+    outfile = os.path.join(tmp_path, 'test.fits')
+    hdu = fits.PrimaryHDU(imdata)
+    hdu.header['SATURATE'] = 35000.0
+    hdu.writeto(outfile, overwrite=True, output_verify='silentfix')
+
+    error_hdu = image_procs.create_error(outfile, maskhdu, rdnoise)
+
+    rms = 0.5 * (
+            np.percentile(imdata[~maskim.astype(bool)], 84.13)
+            - np.percentile(imdata[~maskim.astype(bool)], 15.86)
+    )
+
+    theoretical_error = np.sqrt(imdata + rms**2 + rdnoise)
+
+    np.testing.assert_array_equal(error_hdu.data, theoretical_error)