Merge pull request #102 from LSSTDESC/lsst-tickets/DM-35647

Lsst tickets/dm 35647

.github/workflows/build.yaml

@@ -45,12 +45,24 @@ jobs:
         run: |
           pip install -v -e .
 
+      - name: Run full chain
+        shell: bash -l {0}
+        run: |
+          nosetests tests/run_full_chain.py --all --debug --detailed-errors --verbose --process-restartworker --with-coverage --cover-package=spectractor
+
       - name: Run nosetests
         shell: bash -l {0}
         run: |
-          python setup.py nosetests
+          nosetests tests/run_tests.py --all --debug --detailed-errors --verbose --process-restartworker --with-coverage --cover-package=spectractor
 
-      - name: Run full chain
+      - name: Run doctests and coverage
         shell: bash -l {0}
         run: |
-          nosetests tests/run_full_chain.py --all --debug --detailed-errors --verbose --process-restartworker --with-coverage --cover-package=spectractor
+          ./coverage_and_test.sh
+
+      - name: Coveralls
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+        run: |
+          pip3 install --upgrade coveralls
+          coveralls --service=github

docs/requirements.txt

@@ -5,4 +5,5 @@ pandas
 astropy>3.1
 coloredlogs
 recommonmark
-m2r
+m2r
+mistune==2.0.3

requirements.txt

@@ -22,4 +22,3 @@ deprecated
 pyyaml
 nose
 getCalspec
-pysynphot

spectractor/astrometry.py

@@ -435,9 +435,8 @@ def load_gaia_catalog_around_target(self):
         >>> a = Astrometry("./tests/data/reduc_20170530_134.fits", target_label="HD111980",
         ...                wcs_file_name="./tests/data/reduc_20170530_134_wcs/reduc_20170530_134.wcs")
         >>> a.load_gaia_catalog_around_target()  #doctest: +ELLIPSIS
-        Created TAP+ (v1.2.1) - Connection:...
         INFO: Query finished...
-
+        <Table length=...>...
         """
         radius = 0.5*np.sqrt(2) * max(np.max(self.sources["xcentroid"]) - np.min(self.sources["xcentroid"]),
                                       np.max(self.sources["ycentroid"]) - np.min(self.sources["ycentroid"]))
@@ -1015,7 +1014,7 @@ def run_simple_astrometry(self, extent=None, sources=None):
             # remove background
             self.my_logger.info('\n\tRemove background using astropy SExtractorBackground()...')
             data_wo_bkg = remove_image_background_sextractor(data, sigma=3.0, box_size=(50, 50),
-                                                             filter_size=(10, 10), positive=True)
+                                                             filter_size=(11, 11), positive=True)
             # extract source positions and fluxes
             self.my_logger.info('\n\tDetect sources using photutils source_detection()...')
             self.sources = source_detection(data_wo_bkg, sigma=3.0, fwhm=3.0, threshold_std_factor=5, mask=None)
@@ -1034,7 +1033,8 @@ def run_simple_astrometry(self, extent=None, sources=None):
                   f"--ra {self.target.radec_position.ra.value} --dec {self.target.radec_position.dec.value} " \
                   f"--radius {parameters.CCD_IMSIZE * parameters.CCD_PIXEL2ARCSEC / 3600.} " \
                   f"--dir {self.output_directory} --out {self.tag} " \
-                  f"--overwrite --x-column X --y-column Y {self.sources_file_name}"
+                  f"--overwrite --x-column X --y-column Y {self.sources_file_name} " \
+                  f"--width {self.data.shape[1]} --height {self.data.shape[0]}"
         self.my_logger.info(f'\n\tRun astrometry.net solve_field command:\n\t{command}')
         log = subprocess.check_output(command, shell=True)
         log_file = open(self.log_file_name, "w+")

spectractor/config.py

@@ -79,6 +79,7 @@ def load_config(config_filename, rebin=True):
         FileNotFoundError: Config file ./config/unknown_file.ini does not exist.
 
     """
+    my_logger = set_logger(__name__)
     if not os.path.isfile(os.path.join(parameters.CONFIG_DIR, "default.ini")):
         raise FileNotFoundError('Config file default.ini does not exist.')
     # Load the configuration file
@@ -114,7 +115,7 @@ def load_config(config_filename, rebin=True):
         apply_rebinning_to_parameters()
     else:
         parameters.CCD_REBIN = 1
-        print("No rebinning: parameters.REBIN is forced to 1.")
+        my_logger.warning("No rebinning: parameters.REBIN is forced to 1.")
 
     # check consistency
     if parameters.PIXWIDTH_BOXSIZE > parameters.PIXWIDTH_BACKGROUND:

spectractor/extractor/dispersers.py

@@ -220,15 +220,15 @@ def get_N(deltaX, x0, D=parameters.DISTANCE2CCD, wavelength=656, order=1):
 
 
 def neutral_lines(x_center, y_center, theta_tilt):
-    """Return the nuetrla lines of an hologram."""
+    """Return the neutral lines of a hologram."""
     xs = np.linspace(0, parameters.CCD_IMSIZE, 20)
     line1 = np.tan(theta_tilt * np.pi / 180) * (xs - x_center) + y_center
     line2 = np.tan((theta_tilt + 90) * np.pi / 180) * (xs - x_center) + y_center
     return xs, line1, line2
 
 
 def order01_positions(holo_center, N, theta_tilt, theta0=0, lambda_constructor=639e-6, verbose=True):  # pragma: no cover
-    """Return the order 0 and order 1 positions of an hologram."""
+    """Return the order 0 and order 1 positions of a hologram."""
     # refraction angle between order 0 and order 1 at construction
     alpha = np.arcsin(N * lambda_constructor + np.sin(theta0))
     # distance between order 0 and order 1 in pixels
@@ -250,7 +250,7 @@ def order01_positions(holo_center, N, theta_tilt, theta0=0, lambda_constructor=6
 
 
 def find_order01_positions(holo_center, N_interp, theta_interp, lambda_constructor=639e-6, verbose=True):  # pragma: no cover
-    """Find the order 0 and order 1 positions of an hologram."""
+    """Find the order 0 and order 1 positions of a hologram."""
     N = N_interp(holo_center)
     theta_tilt = theta_interp(holo_center)
     theta0 = 0
@@ -270,13 +270,13 @@ def find_order01_positions(holo_center, N_interp, theta_interp, lambda_construct
 class Grating:
     """Generic class for dispersers."""
 
-    def __init__(self, N, label="", D=parameters.DISTANCE2CCD, data_dir=parameters.DISPERSER_DIR, verbose=False):
+    def __init__(self, N=-1, label="", D=parameters.DISTANCE2CCD, data_dir=parameters.DISPERSER_DIR, verbose=False):
         """Initialize a standard grating object.
 
         Parameters
         ----------
         N: float
-            The number of grooves per mm of the grating
+            The number of grooves per mm of the grating (default: -1)
         label: str
             String label for the grating (default: '')
         D: float
@@ -298,6 +298,8 @@ def __init__(self, N, label="", D=parameters.DISTANCE2CCD, data_dir=parameters.D
         >>> assert g.D is parameters.DISTANCE2CCD
         """
         self.my_logger = set_logger(self.__class__.__name__)
+        if N <= 0 and label == '':
+            raise ValueError("Set either N grooves per mm or the grating label.")
         self.N_input = N
         self.N_err = 1
         self.D = D
@@ -333,7 +335,7 @@ def N(self, x):
         return self.N_input
 
     def load_files(self, verbose=False):
-        """If they exists, load the files in data_dir/label/ to set the main
+        """If they exist, load the files in data_dir/label/ to set the main
         characteristics of the grating. Overrides the N input at initialisation.
 
         Parameters
@@ -366,16 +368,16 @@ def load_files(self, verbose=False):
             a = np.loadtxt(filename)
             self.N_input = a[0]
             self.N_err = a[1]
-        else:
-            raise FileNotFoundError(f"Failed to load {filename} for {self.label}")
+        # else:
+        #     raise FileNotFoundError(f"Failed to load {filename} for {self.label}")
 
         filename = os.path.join(self.data_dir, self.label, "full_name.txt")
         if os.path.isfile(filename):
             with open(filename, 'r') as f:
                 for line in f:  # MFL: you really just want the last line of the file?
                     self.full_name = line.rstrip('\n')
-        else:
-            raise FileNotFoundError(f"Failed to load {filename} for {self.label}")
+        # else:
+        #     raise FileNotFoundError(f"Failed to load {filename} for {self.label}")
 
         filename = os.path.join(self.data_dir, self.label, "transmission.txt")
         if os.path.isfile(filename):
@@ -600,7 +602,7 @@ class Hologram(Grating):
 
     def __init__(self, label, D=parameters.DISTANCE2CCD, data_dir=parameters.DISPERSER_DIR,
                  lambda_plot=256000, verbose=False):
-        """Initialize an Hologram object, given its label. Specification are loaded from text files
+        """Initialize a Hologram object, given its label. Specification are loaded from text files
         in data_dir/label/... Inherit from the Grating class.
 
         Parameters
@@ -734,7 +736,7 @@ def load_specs(self, verbose=True):
         >>> h = Hologram(label='XXX')  # doctest: +ELLIPSIS
         Traceback (most recent call last):
         ...
-        FileNotFoundError:...
+        ValueError:...
 
         """
         if verbose:

spectractor/extractor/images.py

@@ -754,6 +754,7 @@ def find_target(image, guess=None, rotated=False, widths=[parameters.XWINDOW, pa
 
     Examples
     --------
+    >>> parameters.CCD_REBIN = 1
     >>> im = Image('tests/data/reduc_20170605_028.fits', target_label="PNG321.0+3.9")
     >>> im.plot_image()
     >>> guess = [820, 580]

spectractor/extractor/psf.py

@@ -704,7 +704,11 @@ def evaluate(self, pixels, p=None):
         elif pixels.ndim == 1:
             y = np.array(pixels)
             norm = gamma * np.sqrt(np.pi) * special.gamma(alpha - 0.5) / special.gamma(alpha)
-            out = evaluate_moffat1d_unnormalized(y, amplitude, y_c, gamma, alpha) / norm
+            if norm > 0:
+                out = evaluate_moffat1d_unnormalized(y, amplitude, y_c, gamma, alpha)
+                out /= norm
+            else:
+                out = np.zeros_like(y)
             if self.clip:
                 out = np.clip(out, 0, saturation)
             return out

spectractor/extractor/spectrum.py

@@ -522,7 +522,7 @@ def save_spectrum(self, output_file_name, overwrite=False):
             if isinstance(value, astropy.coordinates.angles.Angle):
                 value = value.degree
             hdu1.header[header_key] = value
-            print(f"Set header key {header_key} to {value} from attr {attribute}")
+            # print(f"Set header key {header_key} to {value} from attr {attribute}")
 
         hdu1.header["EXTNAME"] = "SPECTRUM"
         hdu2 = fits.ImageHDU()
@@ -628,7 +628,7 @@ def load_spectrum(self, input_file_name, spectrogram_file_name_override=None,
             for attribute, header_key in fits_mappings.items():
                 if (item := self.header.get(header_key)) is not None:
                     setattr(self, attribute, item)
-                    print(f'set {attribute} to {item}')
+                    # print(f'set {attribute} to {item}')
                 else:
                     print(f'Failed to set spectrum attribute {attribute} using header {header_key}')
 

spectractor/extractor/targets.py

@@ -34,8 +34,8 @@ def load_target(label, verbose=False):
     >>> t = load_target("HD111980", verbose=False)
     >>> print(t.label)
     HD111980
-    >>> print(t.radec_position.dec)
-    -18d31m20.009s
+    >>> print(t.radec_position.dec)  # doctest: +ELLIPSIS
+    -18d31m...s
     >>> parameters.OBS_OBJECT_TYPE = "MONOCHROMATOR"
     >>> t = load_target("XX", verbose=False)
     >>> print(t.label)
@@ -177,8 +177,8 @@ def __init__(self, label, verbose=False):
         >>> s = Star('3C273')
         >>> print(s.label)
         3C273
-        >>> print(s.radec_position.dec)
-        2d03m08.598s
+        >>> print(s.radec_position.dec)  # doctest: +ELLIPSIS
+        2d03m...s
         >>> print(s.emission_spectrum)
         True
 
@@ -187,14 +187,15 @@ def __init__(self, label, verbose=False):
         >>> s = Star('HD111980')
         >>> print(s.label)
         HD111980
-        >>> print(s.radec_position.dec)
-        -18d31m20.009s
+        >>> print(s.radec_position.dec)  # doctest: +ELLIPSIS
+        -18d31m...s
         >>> print(s.emission_spectrum)
         False
 
         """
         Target.__init__(self, label, verbose=verbose)
         self.my_logger = set_logger(self.__class__.__name__)
+        self.simbad_table = None
         self.load()
 
     def load(self):
@@ -226,23 +227,23 @@ def load(self):
         if getCalspec.is_calspec(self.label):
             calspec = getCalspec.Calspec(self.label)
             astroquery_label = calspec.Astroquery_Name
-        simbad_table = simbadQuerier.query_object(astroquery_label)
+        self.simbad_table = simbadQuerier.query_object(astroquery_label)
 
-        if simbad_table is not None:
+        if self.simbad_table is not None:
             if self.verbose or True:
-                self.my_logger.info(f'\n\tSimbad:\n{simbad_table}')
-            self.radec_position = SkyCoord(simbad_table['RA'][0] + ' ' + simbad_table['DEC'][0], unit=(u.hourangle, u.deg))
+                self.my_logger.info(f'\n\tSimbad:\n{self.simbad_table}')
+            self.radec_position = SkyCoord(self.simbad_table['RA'][0] + ' ' + self.simbad_table['DEC'][0], unit=(u.hourangle, u.deg))
         else:
             raise RuntimeError(f"Target {self.label} not found in Simbad")
-        self.get_radec_position_after_pm(simbad_table, date_obs="J2000")
-        if not np.ma.is_masked(simbad_table['Z_VALUE']):
-            self.redshift = float(simbad_table['Z_VALUE'])
+        self.get_radec_position_after_pm(date_obs="J2000")
+        if not np.ma.is_masked(self.simbad_table['Z_VALUE']):
+            self.redshift = float(self.simbad_table['Z_VALUE'])
         else:
             self.redshift = 0
         self.load_spectra()
 
     def load_spectra(self):
-        """Load reference spectra from Pysynphot database or NED database.
+        """Load reference spectra from getCalspec database or NED database.
 
         If the object redshift is >0.2, the LAMBDA_MIN and LAMBDA_MAX parameters
         are redshifted accordingly.
@@ -254,7 +255,7 @@ def load_spectra(self):
         [0.0000000e+00 2.5048577e-14 2.4238061e-14 2.4088789e-14]
         >>> s = Star('HD111980')
         >>> print(s.spectra[0][:4])
-        [2.16890002e-13 2.66480010e-13 2.03540011e-13 2.38780004e-13]
+        [2.3621000e-13 2.1016000e-13 2.1632999e-13 2.4676000e-13]
         >>> s = Star('PKS1510-089')
         >>> print(s.redshift)
         0.36
@@ -338,22 +339,27 @@ def load_spectra(self):
                              f"\n\tEmission spectrum ? {self.emission_spectrum}"
                              f"\n\tLines: {[l.label for l in self.lines.lines]}")
 
-    def get_radec_position_after_pm(self, simbad_table, date_obs):
-        target_pmra = simbad_table[0]['PMRA'] * u.mas / u.yr
-        if np.isnan(target_pmra):
-            target_pmra = 0 * u.mas / u.yr
-        target_pmdec = simbad_table[0]['PMDEC'] * u.mas / u.yr
-        if np.isnan(target_pmdec):
-            target_pmdec = 0 * u.mas / u.yr
-        target_parallax = simbad_table[0]['PLX_VALUE'] * u.mas
-        if target_parallax == 0 * u.mas:
-            target_parallax = 1e-4 * u.mas
-        target_coord = SkyCoord(ra=self.radec_position.ra, dec=self.radec_position.dec,
-                                distance=Distance(parallax=target_parallax),
-                                pm_ra_cosdec=target_pmra, pm_dec=target_pmdec, frame='icrs', equinox="J2000",
-                                obstime="J2000")
-        self.radec_position_after_pm = target_coord.apply_space_motion(new_obstime=Time(date_obs))
-        return self.radec_position_after_pm
+    def get_radec_position_after_pm(self, date_obs):
+        if self.simbad_table is not None:
+            target_pmra = self.simbad_table[0]['PMRA'] * u.mas / u.yr
+            if np.isnan(target_pmra):
+                target_pmra = 0 * u.mas / u.yr
+            target_pmdec = self.simbad_table[0]['PMDEC'] * u.mas / u.yr
+            if np.isnan(target_pmdec):
+                target_pmdec = 0 * u.mas / u.yr
+            target_parallax = self.simbad_table[0]['PLX_VALUE'] * u.mas
+            if target_parallax == 0 * u.mas:
+                target_parallax = 1e-4 * u.mas
+            target_coord = SkyCoord(ra=self.radec_position.ra, dec=self.radec_position.dec,
+                                    distance=Distance(parallax=target_parallax),
+                                    pm_ra_cosdec=target_pmra, pm_dec=target_pmdec, frame='icrs', equinox="J2000",
+                                    obstime="J2000")
+            self.radec_position_after_pm = target_coord.apply_space_motion(new_obstime=Time(date_obs))
+            return self.radec_position_after_pm
+        else:
+            self.my_logger.warning("No Simbad table provided: can't apply proper motion correction. "
+                                   "Return original (RA,DEC) coordinates of the object.")
+            return self.radec_position
 
     def build_sed(self, index=0):
         """Interpolate the database reference spectra and return self.sed as a function of the wavelength.
@@ -368,7 +374,7 @@ def build_sed(self, index=0):
         >>> s = Star('HD111980')
         >>> s.build_sed(index=0)
         >>> s.sed(550)
-        array(1.67605113e-11)
+        array(1.67448019e-11)
         """
         if len(self.spectra) == 0:
             self.sed = interp1d(parameters.LAMBDAS, np.zeros_like(parameters.LAMBDAS), kind='linear', bounds_error=False,