Merge pull request #66 from hover2pi/buggy

Buggy

sedkit/__init__.py

@@ -17,7 +17,7 @@
 __version_commit__ = ''
 _regex_git_hash = re.compile(r'.*\+g(\w+)')
 
-__version__ = '1.1.2'
+__version__ = '1.1.3'
 
 # from pkg_resources import get_distribution, DistributionNotFound
 # try:

sedkit/catalog.py

@@ -42,8 +42,7 @@ def __init__(self, name='SED Catalog', marker='circle', color='blue', verbose=Tr
                      'membership', 'reddening', 'fbol', 'fbol_unc', 'mbol',
                      'mbol_unc', 'Lbol', 'Lbol_unc', 'Lbol_sun',
                      'Lbol_sun_unc', 'Mbol', 'Mbol_unc', 'logg', 'logg_unc',
-                     'mass', 'mass_unc', 'Teff', 'Teff_unc', 'Teff_evo',
-                     'Teff_evo_unc', 'Teff_bb', 'SED']
+                     'mass', 'mass_unc', 'Teff', 'Teff_unc', 'SED']
 
         # A master table of all SED results
         self.results = self.make_results_table(self)
@@ -294,16 +293,21 @@ def from_file(self, filepath, run_methods=['find_2MASS'], delimiter=','):
         # Iterate over table
         for row in data:
 
-            # Make the SED
-            s = SED(row['name'], verbose=False)
-            if 'ra' in row and 'dec' in row:
-                s.sky_coords = row['ra']*q.deg, row['dec']*q.deg
+            try:
+
+                # Make the SED
+                s = SED(row['name'], verbose=False)
+                if 'ra' in row and 'dec' in row:
+                    s.sky_coords = row['ra']*q.deg, row['dec']*q.deg
+
+                # Run the desired methods
+                s.run_methods(run_methods)
 
-            # Run the desired methods
-            s.run_methods(run_methods)
+                # Add it to the catalog
+                self.add_SED(s)
 
-            # Add it to the catalog
-            self.add_SED(s)
+            except:
+                self.message("Could not add SED '{}".format(row['name']))
 
     def get_data(self, *args):
         """Fetch the data for the given columns
@@ -392,13 +396,10 @@ def make_results_table(self):
         results['mass_unc'].unit = q.Msun
         results['Teff'].unit = q.K
         results['Teff_unc'].unit = q.K
-        results['Teff_bb'].unit = q.K
-        results['Teff_evo'].unit = q.K
-        results['Teff_evo_unc'].unit = q.K
 
         return results
 
-    def message(self, msg, pre='[sedkit.Catalog]'):
+    def message(self, msg, pre='[sedkit]'):
         """
         Only print message if verbose=True
 
@@ -416,8 +417,8 @@ def message(self, msg, pre='[sedkit.Catalog]'):
                 print("{} {}".format(pre, msg))
 
     def plot(self, x, y, marker=None, color=None, scale=['linear','linear'],
-             xlabel=None, ylabel=None, fig=None, order=None, identify=None,
-             id_color='red', label_points=False, exclude=None, draw=True, **kwargs):
+             xlabel=None, ylabel=None, fig=None, order=None, identify=[],
+             id_color='red', label_points=False, draw=True, **kwargs):
         """Plot parameter x versus parameter y
 
         Parameters
@@ -458,9 +459,9 @@ def plot(self, x, y, marker=None, color=None, scale=['linear','linear'],
 
         # If no uncertainty column for parameter, add it
         if '{}_unc'.format(x) not in source.column_names:
-            _ = source.add([None] * len(self.source.data['name']), '{}_unc'.format(x))
+            _ = source.add([None] * len(source.data['name']), '{}_unc'.format(x))
         if '{}_unc'.format(y) not in source.column_names:
-            _ = source.add([None] * len(self.source.data['name']), '{}_unc'.format(y))
+            _ = source.add([None] * len(source.data['name']), '{}_unc'.format(y))
 
         # Check if the x parameter is a color
         if '-' in x and all([i in params for i in x.split('-')]):
@@ -500,35 +501,36 @@ def plot(self, x, y, marker=None, color=None, scale=['linear','linear'],
             title = '{} v {}'.format(x, y)
             fig = figure(plot_width=800, plot_height=500, title=title, y_axis_type=scale[1], x_axis_type=scale[0], tools=TOOLS)
 
-        # # Exclude sources
-        # if exclude is not None:
-        #     exc_idx = [i for i, v in enumerate(source.data['name']) if v in exclude]
-        #     patches = {x : [(i, np.nan) for i in exc_idx],
-        #                y : [(i, np.nan) for i in exc_idx]}
-        #
-        #     source.patch(patches)
-
         # Get marker class
         size = kwargs.get('size', 8)
         kwargs['size'] = size
         marker = getattr(fig, marker or self.marker)
         color = color or self.color
-        marker(x, y, source=source, color=color, fill_alpha=0.7, name='points', **kwargs)
 
-        # Plot y errorbars
+        # Prep data
+        names = source.data['name']
+        xval, xerr = source.data[x], source.data['{}_unc'.format(x)]
+        xval[xval == None] = np.nan
+        xerr[xerr == None] = np.nan
         yval, yerr = source.data[y], source.data['{}_unc'.format(y)]
         yval[yval == None] = np.nan
         yerr[yerr == None] = np.nan
+
+        # Plot nominal values
+        marker(x, y, source=source, color=color, fill_alpha=0.7, name='points', **kwargs)
+
+        # Identify sources
+        idx = [ni for ni, name in enumerate(names) if name in identify]
+        fig.circle(xval[idx], yval[idx], size=size + 5, color=id_color, fill_color=None, line_width=2)
+
+        # Plot y errorbars
         y_err_x = [(i, i) for i in source.data[x]]
-        y_err_y = [(i, j) for i, j in zip(yval - yerr, yval + yerr)]
+        y_err_y = [(yval[n] if np.isnan(i - j) else i - j, yval[n] if np.isnan(i + j) else i + j) for n, (i, j) in enumerate(zip(yval, yerr))]
         fig.multi_line(y_err_x, y_err_y, color=color)
 
         # Plot x errorbars
-        xval, xerr = source.data[x], source.data['{}_unc'.format(x)]
-        xval[xval == None] = np.nan
-        xerr[xerr == None] = np.nan
         x_err_y = [(i, i) for i in source.data[y]]
-        x_err_x = [(i, j) for i, j in zip(xval - xerr, xval + xerr)]
+        x_err_x = [(xval[n] if np.isnan(i - j) else i - j, xval[n] if np.isnan(i + j) else i + j) for n, (i, j) in enumerate(zip(xval, xerr))]
         fig.multi_line(x_err_x, x_err_y, color=color)
 
         # Label points
@@ -554,7 +556,7 @@ def plot(self, x, y, marker=None, color=None, scale=['linear','linear'],
 
                 if yerr is not None:
                     ye = np.array(yerr, dtype=float)[idx]
-                    coeffs, cov = np.polyfit(x=xd, y=yd, deg=order, w=1./ye, cov=True)
+                    coeffs, cov = np.polyfit(x=xd, y=yd, deg=order, w=1. / ye, cov=True)
                 else:
                     coeffs, cov = np.polyfit(x=xd, y=yd, deg=order, cov=True)
 
@@ -564,7 +566,7 @@ def plot(self, x, y, marker=None, color=None, scale=['linear','linear'],
                 else:
 
                     # Calculate values and 1-sigma
-                    TT = np.vstack([xaxis**(order-i) for i in range(order + 1)]).T
+                    TT = np.vstack([xaxis**(order - i) for i in range(order + 1)]).T
                     yaxis = np.dot(TT, coeffs)
                     C_yi = np.dot(TT, np.dot(cov, TT.T))
                     sig = np.sqrt(np.diag(C_yi))
@@ -583,24 +585,12 @@ def plot(self, x, y, marker=None, color=None, scale=['linear','linear'],
         # Set axis labels
         xunit = source.data[x].unit
         yunit = source.data[y].unit
-        fig.xaxis.axis_label = '{}{}'.format(x, ' [{}]'.format(xunit) if xunit else '')
-        fig.yaxis.axis_label = '{}{}'.format(y, ' [{}]'.format(yunit) if yunit else '')
+        fig.xaxis.axis_label = xlabel or '{}{}'.format(x, ' [{}]'.format(xunit) if xunit else '')
+        fig.yaxis.axis_label = ylabel or '{}{}'.format(y, ' [{}]'.format(yunit) if yunit else '')
 
         # Formatting
         fig.legend.location = "top_right"
 
-        # Identify sources
-        if isinstance(identify, list):
-            id_cat = Catalog('Identified')
-            for obj_id in identify:
-                obj_result = self.get_SED(obj_id)
-                if str(type(obj_result)) != "<class 'astropy.table.column.Column'>":
-                    obj_result = [obj_result]
-                for obj in obj_result:
-                    id_cat.add_SED(obj)
-            fig = id_cat.plot(x, y, fig=fig, size=size+5, marker='circle', line_color=id_color, fill_color=None, line_width=2, label_points=True)
-            del id_cat
-
         if draw:
             show(fig)
 

sedkit/data/dwarf_sequence.txt

@@ -19,7 +19,7 @@
 # BtVt, Cousins RcIc, Gaia DR2 G/Bp/Rp, Sloan izY, 2MASS JHKs,
 # and WISE W1/W2/W3/W4 photometry. G-V color is Gaia DR2 G - Johnson V.
 #
-SpT   Teff  logT   logL  Mbol  BCv    Mv    B-V    Bt-Vt  G-V    Bp-Rp  G-Rp   M_G    b-y    U-B    V-Rc   V-Ic   V-Ks   J-H    H-Ks    Ks-W1   W1-W2  W1-W3  W1-W4   M_J    M_Ks  i-z  z-Y  R_Rsun Msun  #SpT
+SpT   Teff  logT   Lbol  Mbol  BCv    Mv    B-V    Bt-Vt  G-V    Bp-Rp  G-Rp   M_G    b-y    U-B    V-Rc   V-Ic   V-Ks   J-H    H-Ks    Ks-W1   W1-W2  W1-W3  W1-W4   M_J    M_Ks  i-z  z-Y  radius mass  #SpT
 O3V    46000 4.663  5.80 -9.75 -4.05  -5.7  -0.32   ...    ...    ...    ...    ...    ..... -1.22   ...    .....  .....  .....  .....   ...     ...    ...    ...     ....   ....  ...  ...  12.5   ....  O3V
 O4V    43000 4.633  5.67 -9.42 -3.92  -5.5  -0.32   ...    ...    ...    ...    ...    ..... -1.20   ...    .....  .....  .....  .....   ...     ...    ...    ...     ....   ....  ...  ...  12.3   ....  O4V
 O5V    41500 4.618  5.58 -9.20 -3.77  -5.4  -0.32   ...    ...    ...    ...    ...   -0.133 -1.19   ...    .....  .....  .....  .....   ...     ...    ...    ...     ....   ....  ...  ...  11.9   ....  O5V

sedkit/data/sources.txt

@@ -7,10 +7,8 @@ LHS 1047AB
 LHS 3
 LP 944-20
 LHS 2065
-LHS 1070 ABC
 LHS 1075
 LP 465-27
-LHS 114 AB
 LHS 115
 LP 405-39
 LP 406-29

sedkit/isochrone.py

@@ -72,7 +72,7 @@ def __init__(self, name, units=None, verbose=True, **kwargs):
 
         # Calculate radii if not in the table (R = sqrt(GM/g))
         if 'radius' not in self.data.colnames:
-            radius = np.sqrt((ac.G * (self.data['mass'] * q.M_sun)) / ((10**self.data['logg']) * q.cm / q.s**2)).to(q.R_sun)
+            radius = np.sqrt((ac.G * (self.data['mass'] * q.M_sun)) / ((10**self.data['logg']) * q.m / q.s**2)).to(q.R_sun)
             self.data.add_column(radius, name='radius')
 
         # Get the units
@@ -180,7 +180,7 @@ def evaluate(self, xval, age, xparam, yparam, plot=False):
         if upper is None:
             return None
 
-        # Caluclate the symmetric error
+        # Calculate the symmetric error
         error = max(abs(nominal - lower), abs(nominal - upper)) * 2
 
         # Plot the figure and evaluated point

sedkit/modelgrid.py

@@ -502,7 +502,7 @@ def message(self, msg, pre='[sedkit]'):
             else:
                 print("{} {}".format(pre, msg))
 
-    def photometry(self, bandpasses, weight=True):
+    def photometry(self, bandpasses, weight=False):
         """
         Generate a new ModelGrid object of photometry in the given bands
 
@@ -543,16 +543,20 @@ def photometry(self, bandpasses, weight=True):
 
                 # Calculate the synthetic flux
                 flx, flx_unc = spec.synthetic_flux(bp)
+                wav = bp.wave_eff.astype(np.float16).value
                 if flx is not None and not np.isnan(flx):
                     if err:
-                        data.append([bp.wave_eff.to(self.wave_units).astype(np.float16).value, flx.to(self.flux_units).value, flx_unc.to(self.flux_units).value])
+                        # data.append([bp.wave_eff.astype(np.float16).value, flx.to(self.flux_units).value, flx_unc.to(self.flux_units).value])
+                        data.append([wav, flx.value, flx_unc.value])
                     else:
-                        data.append([bp.wave_eff.to(self.wave_units).astype(np.float16).value, flx.to(self.flux_units).value])
+                        # data.append([bp.wave_eff.astype(np.float16).value, flx.to(self.flux_units).value])
+                        data.append([wav, flx.value])
 
             # Unpack and save as new spectrum
             dataT = np.array(data).T
             phot.add_model(dataT if err else dataT[:2], weights=weights, label=row['label'], **params)
             phot.phot = True
+            phot.native_wave_units = bp.wave_units
 
         return phot
 
@@ -702,23 +706,6 @@ def __init__(self, root=None, **kwargs):
         self.load(root)
 
 
-class Filippazzo2016(ModelGrid):
-    """Child class for the Filippazzo et al. (2016) sample"""
-    def __init__(self):
-        """Load the model object"""
-        model_path = 'data/models/atmospheric/Filippazzo2016.p'
-        root = resource_filename('sedkit', model_path)
-
-        data = pickle.load(open(root, 'rb'))
-
-        # Inherit from base class
-        super().__init__(data['name'], data['parameters'], ref='2017yCat..18100158F')
-
-        # Copy to new __dict__
-        for key, val in data.items():
-            setattr(self, key, val)
-
-
 class SpexPrismLibrary(ModelGrid):
     """Child class for the SpeX Prism Library model grid"""
     def __init__(self):