Merge pull request #82 from mdecleir/plot_tweaks

Plot tweaks

measure_extinction/extdata.py

@@ -949,6 +949,7 @@ def read(self, ext_filename):
                 "ASYM10",
                 "AV11",
                 "AV03",
+                "AV07",
             ]
             self.model["type"] = hdr["MOD_TYPE"]
             for paramkey in paramkeys:
@@ -1111,7 +1112,7 @@ def plot(
             plot x axis as 1/wavelength as is standard for UV extinction curves
 
         exclude : list of strings [default=[]]
-            List of data type(s) to exclude from the plot (e.g., IRS)
+            List of data type(s) to exclude from the plot (e.g., "IRS", "IRAC1",...)
 
         normval : float [default=1.0]
             Normalization value
@@ -1170,6 +1171,10 @@ def plot(
                 x = 1.0 / x
 
             if curtype == "BAND":
+                # do not plot the excluded band(s)
+                for i, bandname in enumerate(self.names[curtype]):
+                    if bandname in exclude:
+                        y[i] = np.nan
                 # plot band data as points with errorbars
                 pltax.errorbar(
                     x, y, yerr=yu, fmt="o", color=color, alpha=alpha, mfc="white"

measure_extinction/merge_obsspec.py

@@ -269,13 +269,15 @@ def merge_spex_obsspec(obstable, mask=[], output_resolution=2000):
     npts[obstable["FLAG"] != 0] = 0
     # take out data points with NaN fluxes
     npts[np.isnan(fluxes)] = 0
+    # quadratically add 1 percent uncertainty to account for unknown uncertainties
+    uncs = np.sqrt(uncs ** 2 + (0.01 * fluxes) ** 2)
     # take out data points with low SNR
     npts[np.less(fluxes / uncs, 10, where=~np.isnan(fluxes / uncs))] = 0
     # take out wavelength regions affected by the atmosphere
-    npts[np.logical_and(1.354e4 < waves, waves < 1.411e4)] = 0
-    npts[np.logical_and(1.805e4 < waves, waves < 1.947e4)] = 0
-    npts[np.logical_and(2.522e4 < waves, waves < 2.875e4)] = 0
-    npts[np.logical_and(4.014e4 < waves, waves < 4.594e4)] = 0
+    npts[np.logical_and(1.347e4 < waves, waves < 1.415e4)] = 0
+    npts[np.logical_and(1.798e4 < waves, waves < 1.949e4)] = 0
+    npts[np.logical_and(2.514e4 < waves, waves < 2.880e4)] = 0
+    npts[np.logical_and(4.000e4 < waves, waves < 4.594e4)] = 0
     # take out data points that need to be masked
     for region in mask:
         npts[(waves > region[0] * 1e4) & (waves < region[1] * 1e4)] = 0

measure_extinction/plotting/plot_ext.py

@@ -21,6 +21,7 @@ def plot_average(
     ax=None,
     extmodels=False,
     fitmodel=False,
+    res=False,
     HI_lines=False,
     range=None,
     exclude=[],
@@ -50,14 +51,17 @@ def plot_average(
     fitmodel: boolean [default=False]
         Whether or not to overplot a fitted model
 
+    res : boolean [default=False]
+        Whether or not to plot the residuals of the fitting (only useful when fitmodel=True)
+
     HI_lines : boolean [default=False]
         Whether or not to indicate the HI-lines in the plot
 
     range : list of 2 floats [default=None]
         Wavelength range to be plotted (in micron) - [min,max]
 
     exclude : list of strings [default=[]]
-        List of data type(s) to exclude from the plot (e.g., IRS)
+        List of data type(s) to exclude from the plot (e.g., "IRS", "IRAC1")
 
     log : boolean [default=False]
         Whether or not to plot the wavelengths on a log-scale
@@ -95,8 +99,8 @@ def plot_average(
     # make a new plot if requested
     if pdf:
         # plotting setup for easier to read plots
-        fontsize = 18
-        font = {"size": fontsize}
+        fs = 20
+        font = {"size": fs}
         plt.rc("font", **font)
         plt.rc("lines", linewidth=1)
         plt.rc("axes", linewidth=2)
@@ -105,9 +109,11 @@ def plot_average(
         plt.rc("ytick.major", width=2, size=10)
         plt.rc("ytick.minor", width=1, size=5)
         plt.rc("axes.formatter", min_exponent=2)
+        plt.rc("xtick", direction="in", labelsize=fs * 0.8)
+        plt.rc("ytick", direction="in", labelsize=fs * 0.8)
 
         # create the plot
-        fig, ax = plt.subplots(figsize=(13, 10))
+        fig, ax = plt.subplots(figsize=(10, 7))
         average.plot(ax, exclude=exclude, color="k")
 
         # plot Milky Way extinction models if requested
@@ -116,7 +122,7 @@ def plot_average(
 
         # overplot a fitted model if requested
         if fitmodel:
-            plot_fitmodel(average, res=True)
+            plot_fitmodel(average, res=res)
 
         # plot HI-lines if requested
         if HI_lines:
@@ -129,10 +135,8 @@ def plot_average(
         # finish configuring the plot
         if log:
             ax.set_xscale("log")
-        plt.xlabel(r"$\lambda$ [$\mu m$]", fontsize=1.5 * fontsize)
-        ax.set_ylabel(
-            average._get_ext_ytitle(ytype=average.type), fontsize=1.5 * fontsize
-        )
+        plt.xlabel(r"$\lambda$ [$\mu m$]", fontsize=fs)
+        ax.set_ylabel(average._get_ext_ytitle(ytype=average.type), fontsize=fs)
         fig.savefig(path + "average_ext.pdf", bbox_inches="tight")
 
         # return the figure and axes for additional manipulations
@@ -234,20 +238,22 @@ def plot_fitmodel(extdata, yoffset=0, res=False):
                 extdata.model["params"][3].value,
             )
         elif extdata.model["type"] == "pow_alav":
-            labeltxt = r"$%5.2f \lambda ^{-%5.2f}$" % (
+            labeltxt = r"$%5.2f \,\lambda^{-%5.2f}$" % (
                 extdata.model["params"][0].value,
                 extdata.model["params"][2].value,
             )
+
         else:
             labeltxt = "fitted model"
         plt.plot(
             extdata.model["waves"],
             extdata.model["exts"] + yoffset,
             "-",
-            lw=2,
-            color="firebrick",
+            lw=3,
+            color="crimson",
             alpha=0.8,
             label=labeltxt,
+            zorder=5,
         )
 
         # plot the underlying powerlaw if a Drude was fitted
@@ -267,8 +273,11 @@ def plot_fitmodel(extdata, yoffset=0, res=False):
 
         # plot the residuals if requested
         if res:
+            plt.setp(plt.gca().get_xticklabels(), visible=False)
             plt.axes([0.125, 0, 0.775, 0.11], sharex=plt.gca())
-            plt.scatter(extdata.model["waves"], extdata.model["residuals"], s=0.5)
+            plt.scatter(
+                extdata.model["waves"], extdata.model["residuals"], s=0.5, color="k"
+            )
             plt.axhline(ls="--", c="k", alpha=0.5)
             plt.axhline(y=0.05, ls=":", c="k", alpha=0.5)
             plt.axhline(y=-0.05, ls=":", c="k", alpha=0.5)
@@ -377,6 +386,7 @@ def plot_multi_extinction(
     log=False,
     text_offsets=[],
     text_angles=[],
+    multicolor=False,
     wavenum=False,
     pdf=False,
 ):
@@ -413,7 +423,7 @@ def plot_multi_extinction(
         Whether or not to spread the extinction curves out by adding a vertical offset to each curve
 
     exclude : list of strings [default=[]]
-        List of data type(s) to exclude from the plot (e.g., IRS)
+        List of data type(s) to exclude from the plot (e.g., "IRS", "IRAC1")
 
     log : boolean [default=False]
         Whether or not to plot the wavelengths on a log-scale
@@ -424,6 +434,8 @@ def plot_multi_extinction(
     text_angles : list of integers [default=[]]
         List of the same length as starpair_list with rotation angles for the annotated text
 
+    multicolor : boolean [default=False]
+        Whether or not to give all curves a different color
     wavenum : boolean [default=False]
         Whether or not to plot the wavelengths as wavenumbers = 1/wavelength
 
@@ -470,8 +482,13 @@ def plot_multi_extinction(
             yoffset = 0.0
 
         # determine where to add the name of the star
-        # find the shortest plotted wavelength
-        (waves, exts, ext_uncs) = extdata.get_fitdata(extdata.waves.keys() - exclude)
+        # find the shortest plotted wavelength, and give preference to spectral data when available
+        exclude2 = []
+        if "BAND" in extdata.waves.keys() and len(extdata.waves.keys()) > 1:
+            exclude2 = ["BAND"]
+        (waves, exts, ext_uncs) = extdata.get_fitdata(
+            extdata.waves.keys() - (exclude + exclude2)
+        )
         if range is not None:
             waves = waves[waves.value >= range[0]]
         min_wave = waves[-1]
@@ -485,9 +502,13 @@ def plot_multi_extinction(
         ann_range = [min_wave, min_wave] * u.micron
 
         # plot the extinction curve
+        if multicolor:
+            pcolor = colors(i % 10)
+        else:
+            pcolor = "k"
         extdata.plot(
             ax,
-            color=colors(i % 10),
+            color=pcolor,
             alpha=0.7,
             alax=alax,
             exclude=exclude,
@@ -606,7 +627,7 @@ def plot_extinction(
         Wavelength range to be plotted (in micron) - [min,max]
 
     exclude : list of strings [default=[]]
-        List of data type(s) to exclude from the plot (e.g., IRS)
+        List of data type(s) to exclude from the plot (e.g., "IRS", "IRAC1")
 
     log : boolean [default=False]
         Whether or not to plot the wavelengths on a log scale
@@ -634,7 +655,7 @@ def plot_extinction(
     plt.rc("axes.formatter", min_exponent=2)
 
     # create the plot
-    fig, ax = plt.subplots(figsize=(13, 10))
+    fig, ax = plt.subplots(figsize=(10, 7))
 
     # read in and plot the extinction curve data for this star
     extdata = ExtData("%s%s_ext.fits" % (path, starpair.lower()))
@@ -695,6 +716,9 @@ def plot_extinction(
     else:
         plt.show()
 
+    # return the figure and axes for additional manipulations
+    return fig, ax
+
 
 def main():
     # commandline parser

measure_extinction/plotting/plot_spec.py

@@ -106,6 +106,7 @@ def plot_multi_spectra(
     class_offset=True,
     text_offsets=[],
     text_angles=[],
+    multicolor=False,
     wavenum=False,
     deredden=False,
     pdf=False,
@@ -153,6 +154,9 @@ def plot_multi_spectra(
     text_angles : list of integers [default=[]]
         List of the same length as starlist with rotation angles for the annotated text
 
+    multicolor : boolean [default=False]
+        Whether or not to give all spectra a different color
+
     pdf : boolean [default=False]
         Whether or not to save the figure as a pdf file
 
@@ -239,9 +243,13 @@ def plot_multi_spectra(
         ann_range = [min_x, min_x] * u.micron
 
         # plot the spectrum
+        if multicolor:
+            pcolor = colors(i % 10)
+        else:
+            pcolor = "k"
         starobs.plot(
             ax,
-            pcolor=colors(i % 10),
+            pcolor=pcolor,
             norm_wave_range=norm_range,
             mlam4=mlam4,
             wavenum=wavenum,
@@ -253,7 +261,7 @@ def plot_multi_spectra(
             annotate_text=star.upper() + "  " + starobs.sptype,
             annotate_yoffset=text_offsets[i],
             annotate_rotation=text_angles[i],
-            annotate_color=colors(i % 10),
+            annotate_color=pcolor,
         )
 
     # plot HI-lines if requested
@@ -426,6 +434,9 @@ def plot_spectrum(
     else:
         plt.show()
 
+    # return the figure and axes for additional manipulations
+    return fig, ax
+
 
 def main():
     # commandline parser

measure_extinction/stardata.py

@@ -1175,7 +1175,7 @@ def plot(
             plot x axis as 1/wavelength as is standard for UV extinction curves
 
         exclude : list of strings [default=[]]
-            Which data type(s) to exclude from the plot (e.g., IRS)
+            List of data type(s) to exclude from the plot (e.g., "IRS", "IRAC1",...)
 
         yoffset : float [default=None]
             multiplicative or additive offset for the data
@@ -1306,6 +1306,10 @@ def plot(
             else:
                 yplotvals += yoffset
             if curtype == "BAND":
+                # do not plot the excluded band(s)
+                for i, bandname in enumerate(self.data[curtype].get_band_names()):
+                    if bandname in exclude:
+                        yplotvals[i] = np.nan
                 # plot band data as points with errorbars
                 ax.errorbar(
                     x,

measure_extinction/utils/merge_spex_spec.py

@@ -5,7 +5,6 @@
 import argparse
 import pkg_resources
 import os
-
 from astropy.table import Table
 
 from measure_extinction.merge_obsspec import merge_spex_obsspec