Fix compatibility issue between old and new versions of the IGRINS PL…

…P for spec_a0v.fits files.

src/muler/hpf.py

@@ -12,12 +12,14 @@
 import warnings
 import logging
 from muler.echelle import EchelleSpectrum, EchelleSpectrumList
+from muler.utilities import resample_list
 import numpy as np
 import astropy
 from astropy.io import fits
 from astropy import units as u
 from astropy.wcs import WCS, FITSFixedWarning
 from astropy.nddata import StdDevUncertainty
+from scipy.ndimage import median_filter
 from scipy.interpolate import InterpolatedUnivariateSpline
 from astropy.constants import R_jup, R_sun, G, M_jup, R_earth, c
 from astropy.time import Time
@@ -28,6 +30,7 @@
 from . import templates
 import pandas as pd
 
+
 log = logging.getLogger(__name__)
 
 for category in [
@@ -349,6 +352,9 @@ def sky_subtract(self, method="scalar", scale=0.93):
             beta_native_spectrum = spec.get_static_sky_ratio_template()
             resampler = LinearInterpolatedResampler(extrapolation_treatment="zero_fill")
             beta = resampler(beta_native_spectrum, spec.spectral_axis)
+        # elif method == 'median': #Attempt to measure the scale of the sky lines between the fibers using a median of the pixels
+        #     resampled_sky = spec.sky.resample(spec)
+        #     good_sci_pixels = self.flux.value -
         else:
             log.error("Method must be one of 'naive', 'scalar' or 'vector'. ")
             raise NotImplementedError
@@ -475,6 +481,24 @@ def sky_subtract(self, method="vector", scale=0.93):
 
         return spec_out
 
+    def test_print_sky_scale(self):
+        #reampled_sky = resample_list(self.sky, self) #Resample sky to match sci wavelengths
+        all_sky = np.zeros([2048, len(self)]) #Put all science and sky fiber flux into 2D arrays for easy manipulation
+        all_sci = np.zeros([2048, len(self)])
+        for i in range(len(self)):
+            all_sky[:, i] = (self[i].sky).resample(self[i]).flux.value
+            all_sci[:, i] = self[i].flux.value
+        all_sky -= median_filter(all_sky, [25, 1]) #Subtract continuum/background using a running median filter
+        all_sci -= median_filter(all_sci, [25, 1])
+        #all_sky[all_sky ]
+        max_sky_flux = np.nanmax(all_sky)
+        bad_pix = (all_sky < 0.1 * max_sky_flux) & (all_sci < 0.1 * max_sky_flux)
+        all_sky[bad_pix] = np.nan
+        all_sci[bad_pix] = np.nan
+        print(np.nanmedian(all_sci / all_sky))
+
+
+
     # def sky_subtract(self):
     #     """Sky subtract all orders
     #     """

src/muler/igrins.py

@@ -494,11 +494,14 @@ def read(file, precache_hdus=True, wavefile=None):
                 wave_hdus = fits.open(full_path, memmap=False)
             cached_hdus.append(wave_hdus)
 
-        hdus0_shape = hdus[0].data.shape
+        if hdus[0].data is not None:
+            hdus0_shape = hdus[0].data.shape #Normally we read from the 0th extension
+        else:
+            hdus0_shape = hdus[1].data.shape #To insure compatibility with new version of the PLP for spec_a0v.fits files
         if len(hdus0_shape) == 2: #1D spectrum
-            n_orders, n_pix = hdus[0].data.shape
+            n_orders, n_pix = hdus0_shape
         elif len(hdus0_shape) == 3: #2D spectrum
-            n_orders, n_height, n_pix = hdus[0].data.shape
+            n_orders, n_height, n_pix = hdus0_shape
 
         list_out = []
         for i in range(n_orders - 1, -1, -1):