Merge pull request #612 from IainHammond/master

new fits.open_header function, typo fixes and better verbosity

vip_hci/fits/fits.py

@@ -4,15 +4,19 @@
 """
 
 
-__author__ = "C. A. Gomez Gonzalez, T. Bédrine, V. Christiaens"
-__all__ = ["open_fits", "info_fits", "write_fits", "verify_fits", 
+__author__ = "C. A. Gomez Gonzalez, T. Bédrine, V. Christiaens, I. Hammond"
+__all__ = ["open_fits", "info_fits", "write_fits", "verify_fits",
            "byteswap_array"]
 
 
-import os
+from os.path import isfile, exists
+from os import remove
+
 import numpy as np
-from astropy.io import fits as ap_fits
-from astropy.io.fits import HDUList
+from astropy.io.fits.convenience import writeto
+from astropy.io.fits.hdu.hdulist import fitsopen, HDUList
+from astropy.io.fits.hdu.image import ImageHDU
+
 from ..config.paramenum import ALL_FITS
 
 
@@ -57,21 +61,19 @@ def open_fits(fitsfilename, n=0, header=False, ignore_missing_end=False,
         equals -2, returns a list of all arrays.
     header : dict or list of dict
         [memmap=False, header=True] Dictionary containing the fits header.
-        If n equals -2, returns a list of all dictionnaries.
+        If n equals -2, returns a list of all dictionaries.
 
     """
     fitsfilename = str(fitsfilename)
-    if not os.path.isfile(fitsfilename):
+    if not isfile(fitsfilename):
         fitsfilename += ".fits"
-    
+
     try:
-        hdulist = ap_fits.open(fitsfilename, 
-                               ignore_missing_end=ignore_missing_end,
-                   	           	memmap=True, **kwargs)
-    except: # If BZERO/BSCALE/BLANK header keywords are present, HDU can’t be loaded as memory map
-        hdulist = ap_fits.open(fitsfilename,
-                               ignore_missing_end=ignore_missing_end,
-                   	           	memmap=False, **kwargs)
+        hdulist = fitsopen(fitsfilename, ignore_missing_end=ignore_missing_end,
+                           memmap=True, **kwargs)
+    except ValueError:  # If BZERO/BSCALE/BLANK header keywords are present, HDU can’t be loaded as memory map
+        hdulist = fitsopen(fitsfilename, ignore_missing_end=ignore_missing_end,
+                           memmap=False, **kwargs)
 
     # Opening all extensions in a MEF
     if n == ALL_FITS:
@@ -82,26 +84,27 @@ def open_fits(fitsfilename, n=0, header=False, ignore_missing_end=False,
 
         for index, element in enumerate(hdulist):
             data, head = _return_data_fits(hdulist=hdulist, index=index,
-                                           precision=precision,
+                                           header=header, precision=precision,
                                            verbose=verbose)
             data_list.append(data)
             header_list.append(head)
 
         hdulist.close()
         if header:
             if verbose:
-                print(f"All fits HDU data and headers succesfully loaded.")
+                print(f"All {len(hdulist)} FITS HDU data and headers successfully loaded. ")
             return data_list, header_list
         else:
             if verbose:
-                print(f"All fits HDU data succesfully loaded.")
+                print(f"All {len(hdulist)} FITS HDU data successfully loaded. ")
             return data_list
+
     # Opening only a specified extension
     else:
         if return_memmap:
             return hdulist[n]
 
-        data, head = _return_data_fits(hdulist=hdulist, index=n, 
+        data, head = _return_data_fits(hdulist=hdulist, index=n, header=header,
                                        precision=precision, verbose=verbose)
         hdulist.close()
         if header:
@@ -112,9 +115,9 @@ def open_fits(fitsfilename, n=0, header=False, ignore_missing_end=False,
 
 def _return_data_fits(hdulist: HDUList,
                       index: int,
+                      header: bool = False,
                       precision=np.float32,
-                      verbose: bool = True,
-):
+                      verbose: bool = True):
     """
     Subfunction used to return data (and header) from a given index.
 
@@ -130,10 +133,13 @@ def _return_data_fits(hdulist: HDUList,
     head = hdulist[index].header
 
     if verbose:
-        print(
-            f"Fits HDU-{index} data successfully loaded, header available. "
-            f"Data shape: {data.shape}"
-        )
+        if header:
+            print(f"FITS HDU-{index} data and header successfully loaded. "
+                  f"Data shape: {data.shape}")
+        else:
+            print(f"FITS HDU-{index} data successfully loaded. "
+                  f"Data shape: {data.shape}")
+
     return data, head
 
 
@@ -163,7 +169,7 @@ def byteswap_array(array):
     Note
     ----
     More info about byteswapping here:
-    http://docs.scipy.org/doc/numpy-1.10.1/user/basics.byteswapping.html
+    https://docs.scipy.org/doc/numpy-1.10.1/user/basics.byteswapping.html
 
     """
     array_out = array.byteswap().newbyteorder()
@@ -183,7 +189,7 @@ def info_fits(fitsfilename, **kwargs):
         "output_verify" can be set to ignore, in case of non-standard header.
 
     """
-    with ap_fits.open(fitsfilename, memmap=True, **kwargs) as hdulist:
+    with fitsopen(fitsfilename, memmap=True, **kwargs) as hdulist:
         hdulist.info()
 
 
@@ -199,10 +205,10 @@ def verify_fits(fitsfilename):
     """
     if isinstance(fitsfilename, list):
         for ffile in fitsfilename:
-            with ap_fits.open(ffile) as f:
+            with fitsopen(ffile) as f:
                 f.verify()
     else:
-        with ap_fits.open(fitsfilename) as f:
+        with fitsopen(fitsfilename) as f:
             f.verify()
 
 
@@ -238,12 +244,12 @@ def write_fits(fitsfilename, array, header=None, output_verify="exception",
         fitsfilename += ".fits"
 
     res = "saved"
-    if os.path.exists(fitsfilename):
-        os.remove(fitsfilename)
+    if exists(fitsfilename):
+        remove(fitsfilename)
         res = "overwritten"
 
     if isinstance(array, tuple):
-        new_hdul = ap_fits.HDUList()
+        new_hdul = HDUList()
         if header is None:
             header = [None] * len(array)
         elif not isinstance(header, tuple):
@@ -255,12 +261,12 @@ def write_fits(fitsfilename, array, header=None, output_verify="exception",
 
         for i in range(len(array)):
             array_tmp = array[i].astype(precision, copy=False)
-            new_hdul.append(ap_fits.ImageHDU(array_tmp, header=header[i]))
+            new_hdul.append(ImageHDU(array_tmp, header=header[i]))
 
         new_hdul.writeto(fitsfilename, output_verify=output_verify)
     else:
         array = array.astype(precision, copy=False)
-        ap_fits.writeto(fitsfilename, array, header, output_verify)
+        writeto(fitsfilename, array, header, output_verify)
 
     if verbose:
-        print("Fits file successfully {}".format(res))
+        print(f"Fits file successfully {res}")

vip_hci/fits/headers.py

@@ -3,13 +3,16 @@
 Module with conversion utilities from dictionaries to headers, and reversely.
 """
 
-__author__ = "Thomas Bedrine"
+__author__ = "Thomas Bedrine, Iain Hammond"
 __all__ = ["dict_to_fitsheader",
-           "fitsheader_to_dict"]
-
+           "fitsheader_to_dict",
+           "open_header"]
 
+from os.path import isfile
 from typing import Tuple
-from astropy.io.fits import Header
+
+from astropy.io.fits.convenience import getheader
+from astropy.io.fits.header import Header
 
 
 def dict_to_fitsheader(initial_dict: dict) -> Header:
@@ -37,7 +40,7 @@ def fitsheader_to_dict(
     initial_header: Header, sort_by_prefix: str = ""
 ) -> Tuple[dict, str]:
     """
-    Extract a dictionnary of parameters and a string from a FITS Header.
+    Extract a dictionary of parameters and a string from a FITS Header.
 
     The string is supposedly the name of the algorithm that was used to obtain
     the results that go with the Header.
@@ -71,3 +74,48 @@ def fitsheader_to_dict(
     algo_name = parameters["algo_name"]
     del parameters["algo_name"]
     return parameters, algo_name
+
+
+def open_header(fitsfilename: str, n: int = 0, extname: str = None,
+                verbose: bool = False) -> Header:
+    """
+    Load a FITS header into memory to avoid loading the data.
+
+    This function is a simple wrapper of astropy.io.fits.convenience.getheader
+    designed to substitute `open_fits` when only a FITS header is needed.
+    This is ~ 40 times faster than using `open_fits` with header=True on
+    an average sized VLT/SPHERE data set.
+
+    Parameters
+    ----------
+    fitsfilename : string
+        Name of the FITS file.
+    n : int, optional
+        Which HDU ext to open. Default is the first ext (zero based indexing).
+    extname : str, optional
+        Opens the HDU ext by name, rather than by HDU number. Overrides `n` and
+        is not case-sensitive.
+    verbose : bool, optional
+        If True prints message of completion.
+
+    Returns
+    -------
+    header : dict or list of dict
+        Astropy Header class with both a dict-like and list-like interface.
+    """
+
+    fitsfilename = str(fitsfilename)
+    if not isfile(fitsfilename):
+        fitsfilename += ".fits"
+
+    # if extname is a non-empty string, provide that instead. Else use HDU number
+    if extname and not extname.isspace():
+        n = extname
+        header = getheader(fitsfilename, extname=n, ignore_missing_end=True)
+    else:
+        header = getheader(fitsfilename, ext=n, ignore_missing_end=True)
+
+    if verbose:
+        print(f"Fits HDU-{n} header successfully loaded.")
+
+    return header

vip_hci/fm/negfc_fmerit.py

@@ -68,7 +68,7 @@ def chisquare(
     psfs_norm: numpy.array
         The scaled psf expressed as a numpy.array.
     fwhm : float
-        The FHWM in pixels.
+        The FWHM in pixels.
     annulus_width: int, optional
         The width in pixels of the annulus on which the PCA is done.
     aperture_radius: int, optional

vip_hci/fm/negfc_mcmc.py

@@ -137,7 +137,7 @@ def lnlike(param, cube, angs, psf_norm, fwhm, annulus_width, ncomp,
     ncomp: int or None
         The number of principal components for PCA-based algorithms.
     fwhm : float
-        The FHWM in pixels.
+        The FWHM in pixels.
     aperture_radius: float
         The radius of the circular aperture in terms of the FWHM.
     initial_state: numpy.array
@@ -335,7 +335,7 @@ def lnprob(param, bounds, cube, angs, psf_norm, fwhm, annulus_width, ncomp,
         If the input cube is 4D, psfn must be either 3D or 4D. In either cases,
         the first dimension(s) must match those of the input cube.
     fwhm : float
-        The FHWM in pixels.
+        The FWHM in pixels.
     annulus_width: float
         The width in pixel of the annulus on wich the PCA is performed.
     ncomp: int or None
@@ -525,7 +525,7 @@ def mcmc_negfc_sampling(cube, angs, psfn, initial_state, algo=pca_annulus,
     aperture_radius: float, optional
         The radius in FWHM of the circular aperture.
     fwhm : float
-        The FHWM in pixels.
+        The FWHM in pixels.
     mu_sigma: tuple of 2 floats or bool, opt
         If set to None: not used, and falls back to original version of the
         algorithm, using ``fmerit`` [WER17]_.

vip_hci/fm/negfc_nested.py

@@ -88,7 +88,7 @@ def nested_negfc_sampling(init, cube, angs, psfn, fwhm, mu_sigma=True,
         If the input cube is 4D, psfn must be either 3D or 4D. In either cases,
         the first dimension(s) must match those of the input cube.
     fwhm : float
-        The FHWM in pixels.
+        The FWHM in pixels.
     mu_sigma: tuple of 2 floats or bool, opt
         If set to None: not used, and falls back to original version of the
         algorithm, using fmerit [WER17]_.

vip_hci/fm/negfc_simplex.py

@@ -54,7 +54,7 @@ def firstguess_from_coord(planet, center, cube, angs, psfn, fwhm, annulus_width,
         If the input cube is 4D, psfn must be either 3D or 4D. In either cases,
         the first dimension(s) must match those of the input cube.
     fwhm : float
-        The FHWM in pixels.
+        The FWHM in pixels.
     annulus_width: int, optional
         The width in pixels of the annulus on which the PCA is done.
     aperture_radius: int, optional
@@ -323,7 +323,7 @@ def firstguess_simplex(p, cube, angs, psfn, ncomp, fwhm, annulus_width,
     ncomp: int or None
         The number of principal components.
     fwhm : float
-        The FHWM in pixels.
+        The FWHM in pixels.
     annulus_width: int, optional
         The width in pixels of the annulus on which the PCA is done.
     aperture_radius: int, optional
@@ -479,7 +479,7 @@ def firstguess(cube, angs, psfn, ncomp, planets_xy_coord, fwhm=4,
     plsc: float, optional
         The platescale, in arcsec per pixel.
     fwhm : float, optional
-        The FHWM in pixels.
+        The FWHM in pixels.
     annulus_width: int, optional
         The width in pixels of the annulus on which the PCA is done.
     aperture_radius: int, optional

vip_hci/metrics/detection.py

@@ -13,7 +13,7 @@
 import numpy as np
 import pandas as pn
 from hciplot import plot_frames
-from scipy.ndimage.filters import correlate
+from scipy.ndimage import correlate
 from skimage import feature
 from astropy.stats import sigma_clipped_stats
 from astropy.stats import gaussian_fwhm_to_sigma, gaussian_sigma_to_fwhm

vip_hci/objects/dataset.py

@@ -601,7 +601,7 @@ def __init__(
         self.fwhm = fwhm
         if self.fwhm is not None:
             if self.cube.ndim == 4:
-                check_array(self.fwhm, dim=1, msg="FHWM")
+                check_array(self.fwhm, dim=1, msg="FWHM")
             elif self.cube.ndim == 3:
                 print("FWHM: {}".format(self.fwhm))
         self.px_scale = px_scale

vip_hci/preproc/badpixremoval.py

@@ -15,7 +15,8 @@
 
 
 import warnings
-__author__ = 'V. Christiaens, Carlos Alberto Gomez Gonzalez'
+__author__ = ('V. Christiaens, Carlos Alberto Gomez Gonzalez, '
+              'Srikanth Kompella')
 __all__ = ['frame_fix_badpix_isolated',
            'cube_fix_badpix_isolated',
            'cube_fix_badpix_annuli',
@@ -334,7 +335,8 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, correct_only=False,
                 array_out[i] = res[0]
                 final_bpm[i] = res[1]
         else:
-            print("Processing using ADACS' multiprocessing approach...")
+            if verbose:
+                print("Cleaning frames using ADACS' multiprocessing approach")
             # dummy calling the function to create cached version of the code prior to forking
             if bpm_mask is not None:
                 bpm_mask_dum = bpm_mask[0]
@@ -1110,8 +1112,8 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori,
                                         verbose)
                     array_corr[i], bpix_map_cumul[i] = res
             else:
-                msg = "Cleaning frames using ADACS' multiprocessing approach"
-                print(msg)
+                if verbose:
+                    print("Cleaning frames using ADACS' multiprocessing approach")
                 # creating shared memory buffer space for the image cube.
                 shm_clump = shared_memory.SharedMemory(create=True,
                                                        size=array_corr.nbytes)
@@ -1182,8 +1184,8 @@ def _mp_clump_slow(args):
                                                  neighbor_box, nneig, half_res_y,
                                                  verbose)
             else:
-                msg = "Cleaning frames using ADACS' multiprocessing approach"
-                print(msg)
+                if verbose:
+                    print("Cleaning frames using ADACS' multiprocessing approach")
                 # dummy calling sigma_filter function to create a cached version of the numba function
                 if bpm_mask.ndim == 3:
                     dummy_bpm = bpm_mask[0]

vip_hci/psfsub/framediff.py

@@ -76,7 +76,7 @@ def frame_diff(*all_args: List, **all_kwargs: dict):
     angle_list : numpy ndarray, 1d
         Corresponding parallactic angle for each frame.
     fwhm : float, optional
-        Known size of the FHWM in pixels to be used. Default is 4.
+        Known size of the FWHM in pixels to be used. Default is 4.
     metric : str, optional
         Distance metric to be used ('cityblock', 'cosine', 'euclidean', 'l1',
         'l2', 'manhattan', 'correlation', etc). It uses the scikit-learn