Skip to content

Commit

Permalink
update plotting, ci fixes (#41)
Browse files Browse the repository at this point in the history 
Nice job!
  • Loading branch information
@alexmaragko
alexmaragko authored Feb 29, 2024
1 parent 3082e96 commit ebcdbba
Show file tree
Hide file tree
Showing 4 changed files with 49 additions and 41 deletions.
59 changes: 33 additions & 26 deletions pypahdb/decomposer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ def save_pdf(self, filename, header="", domaps=True, doplots=True):
d = pdf.infodict()
d["Title"] = "pyPAHdb Results Summary"
d["Author"] = (
"Dr. C. Boersma, Dr. M.J. Shannon, and Dr. A. " "Maragkoudakis"
"Dr. C. Boersma, Dr. M.J. Shannon, and Dr. A. Maragkoudakis"
)
d["Producer"] = "NASA Ames Research Center"
d["Creator"] = "pypahdb v{}(Python {}.{}.{})".format(
Expand Down Expand Up @@ -158,7 +158,7 @@ def _fits_to_disk(hdr, filename):
"Software used to create this file",
)
hdr["AUTHOR"] = (
"Dr. C. Boersma, Dr. M.J. Shannon, and Dr. A. " "Maragkoudakis",
"Dr. C. Boersma, Dr. M.J. Shannon, and Dr. A. Maragkoudakis",
"Authors of the software",
)
cards = [
Expand Down Expand Up @@ -382,11 +382,6 @@ def plot_fit(self, i=0, j=0):
"""

# Enable quantity_support.
from astropy.visualization import quantity_support

quantity_support()

# Create figure on shared axes.
fig = plt.figure()
gs = gridspec.GridSpec(4, 1, height_ratios=[2, 1, 2, 2], figure=fig)
Expand All @@ -397,10 +392,18 @@ def plot_fit(self, i=0, j=0):
ax1 = fig.add_subplot(gs[1], sharex=ax0)
ax2 = fig.add_subplot(gs[2], sharex=ax0)
ax3 = fig.add_subplot(gs[3], sharex=ax0)

# Convenience defintions.
for ax in [ax0, ax2, ax3]:
fmt = ax.yaxis.get_major_formatter()
fmt.set_scientific(True)
fmt.set_powerlimits((-3, 1))
for ax in [ax0, ax1, ax2]:
plt.setp(ax.get_xticklabels(), visible=False)

# Convenience definitions.
abscissa = self.spectrum.spectral_axis
charge = self.charge
# Check if size of datapoints are too large and change marker size.
ms = 5 if len(abscissa) < 1000 else 2

# ax0: Best fit.
data = self.spectrum.flux.T[:, i, j]
Expand All @@ -412,8 +415,8 @@ def plot_fit(self, i=0, j=0):
abscissa,
data,
yerr=unc,
marker="x",
ms=5,
marker=".",
ms=ms,
mew=0.5,
lw=0,
color="black",
Expand All @@ -422,45 +425,47 @@ def plot_fit(self, i=0, j=0):
label="input",
zorder=0,
)
ax0.plot(abscissa, model, label="fit", color="red", lw=1.5)
ax0.plot(abscissa, model, label="fit", color="tab:red", lw=1.5)
error_str = "$error$=%-4.2f" % (self.error[i][j])
ax0.text(0.025, 0.9, error_str, ha="left", va="center", transform=ax0.transAxes)
ax0.text(0.025, 0.88, error_str, ha="left", va="center", transform=ax0.transAxes)
ax0.set_ylabel(f'{self.spectrum.meta["colnames"][1]} [{self.spectrum.flux.unit}]')

# ax1: Residual.
ax1.plot(abscissa, data - model, lw=1, label="residual", color="black")
ax1.plot(abscissa, data - model, lw=1, label="residual", color="gray")
ax1.axhline(y=0, color="0.5", ls="--", dashes=(12, 16), zorder=-10, lw=0.5)

# ax2: Size breakdown.
ax2.errorbar(
abscissa,
data,
yerr=unc,
marker="x",
ms=5,
marker=".",
ms=ms,
mew=0.5,
lw=0,
color="black",
ecolor="grey",
capsize=2,
zorder=0,
)
ax2.plot(abscissa, model, color="red", lw=1.5)
ax2.plot(abscissa, model, color="tab:red", lw=1.5)
ax2.plot(
abscissa, self.size["large"][:, i, j], label="large", lw=1, color="purple"
abscissa, self.size["large"][:, i, j], label="large", lw=1, color="tab:green"
)
ax2.plot(
abscissa, self.size["small"][:, i, j], label="small", lw=1, color="crimson"
abscissa, self.size["small"][:, i, j], label="small", lw=1, color="tab:blue"
)
size_str = "$f_{large}$=%3.1f" % (self.large_fraction[i][j])
ax2.text(0.025, 0.9, size_str, ha="left", va="center", transform=ax2.transAxes)
ax2.text(0.025, 0.88, size_str, ha="left", va="center", transform=ax2.transAxes)
ax2.set_ylabel(f'{self.spectrum.meta["colnames"][1]} [{self.spectrum.flux.unit}]')

# ax3: Charge breakdown.
ax3.errorbar(
abscissa,
data,
yerr=unc,
marker="x",
ms=5,
marker=".",
ms=ms,
mew=0.5,
lw=0,
color="black",
Expand All @@ -470,16 +475,18 @@ def plot_fit(self, i=0, j=0):
)
ax3.plot(abscissa, model, color="red", lw=1.5)
ax3.plot(
abscissa, charge["anion"][:, i, j], label="anion", lw=1, color="orange"
abscissa, charge["anion"][:, i, j], label="anion", lw=1, color="tab:orange"
)
ax3.plot(
abscissa, charge["neutral"][:, i, j], label="neutral", lw=1, color="green"
abscissa, charge["neutral"][:, i, j], label="neutral", lw=1, color="tab:cyan"
)
ax3.plot(
abscissa, charge["cation"][:, i, j], label="cation", lw=1, color="blue"
abscissa, charge["cation"][:, i, j], label="cation", lw=1, color="tab:purple"
)
ion_str = "$f_{ionized}$=%3.1f" % (self.ionized_fraction[i][j])
ax3.text(0.025, 0.9, ion_str, ha="left", va="center", transform=ax3.transAxes)
ax3.text(0.025, 0.88, ion_str, ha="left", va="center", transform=ax3.transAxes)
ax3.set_xlabel(f'{self.spectrum.meta["colnames"][0]} [{self.spectrum.spectral_axis.unit}]')
ax3.set_ylabel(f'{self.spectrum.meta["colnames"][1]} [{self.spectrum.flux.unit}]')

# Set tick parameters and add legends to axes.
for ax in (ax0, ax1, ax2, ax3):
Expand Down
17 changes: 9 additions & 8 deletions pypahdb/observation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ def __init__(self, file_path):
warnings.simplefilter("ignore", category=VerifyWarning)

self.spectrum = Spectrum1D.read(self.file_path)
self.spectrum.meta['colnames'] = 3 * ["",]

# Always work as if spectrum is a cube.
if len(self.spectrum.flux.shape) == 1:
Expand Down Expand Up @@ -92,6 +93,7 @@ def __init__(self, file_path):
flux = h.data.T * u.Unit(h.header["BUNIT"])
wave = hdu[h0].data[h1] * u.Unit(hdu[h0].columns[h1].unit)
self.spectrum = Spectrum1D(flux, spectral_axis=wave)
self.spectrum.meta['colnames'] = 3 * ["",]

return None

Expand Down Expand Up @@ -125,31 +127,30 @@ def __init__(self, file_path):
data = ascii.read(self.file_path)
# Always work as if spectrum is a cube.
flux = np.reshape(
data["FLUX"].quantity,
data[data.colnames[1]].quantity,
(
1,
1,
)
+ data["FLUX"].quantity.shape,
+ data[data.colnames[1]].quantity.shape,
)
# Create Spectrum1D object.
for name in data.colnames:
data.rename_column(name, name.upper())
wave = data["WAVELENGTH"].quantity
wave = data[data.colnames[0]].quantity
unc = None
if "FLUX_UNCERTAINTY" in data.colnames:
if len(data.colnames) > 2:
unc = StdDevUncertainty(
np.reshape(
data["FLUX_UNCERTAINTY"].quantity,
data[data.colnames[2]].quantity,
(
1,
1,
)
+ (data["FLUX_UNCERTAINTY"].quantity.shape),
+ (data[data.colnames[2]].quantity.shape),
)
)

self.spectrum = Spectrum1D(flux, spectral_axis=wave, uncertainty=unc)
self.spectrum.meta['colnames'] = data.colnames

hdr = ""
for card in data.meta["keywords"].keys():
Expand Down
4 changes: 2 additions & 2 deletions pypahdb/tests/test_decomposer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,10 +20,10 @@ class DecomposerTestCase(unittest.TestCase):
"""Unit tests for `decomposer.py`."""

def setUp(self):
import pkg_resources
import importlib_resources
import tempfile
file_name = 'resources/sample_data_NGC7023.tbl'
file_path = pkg_resources.resource_filename('pypahdb', file_name)
file_path = importlib_resources.files('pypahdb') / file_name
self.observation = Observation(file_path)
self.decomposer = Decomposer(self.observation.spectrum)
self.tmpdir = tempfile.gettempdir()
Expand Down
10 changes: 5 additions & 5 deletions pypahdb/tests/test_observation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
"""

import unittest
import pkg_resources
import importlib_resources

from pypahdb.observation import Observation

Expand All @@ -17,21 +17,21 @@ class SpectrumTestCase(unittest.TestCase):
def test_read_spectrum1d(self):
"""Can we create an instance of Observation from a Spectrum1D file?"""
file_name = 'resources/sample_data_jwst.fits'
file_path = pkg_resources.resource_filename('pypahdb', file_name)
file_path = importlib_resources.files('pypahdb') / file_name

assert isinstance(Observation(file_path), Observation)

def test_read_ascii(self):
"""Can we create an instance of Observation from an ASCII file?"""
file_name = 'resources/sample_data_NGC7023.tbl'
file_path = pkg_resources.resource_filename('pypahdb', file_name)
file_path = importlib_resources.files('pypahdb') / file_name

assert isinstance(Observation(file_path), Observation)

def test_read_fits(self):
"""Can we create an instance of Observation from a FITS file?"""
file_name = 'resources/sample_data_NGC7023.fits'
file_path = pkg_resources.resource_filename('pypahdb', file_name)
file_path = importlib_resources.files('pypahdb') / file_name

assert isinstance(Observation(file_path), Observation)

Expand All @@ -44,7 +44,7 @@ def test_file_not_found(self):
def test_file_malformed(self):
"""Can we detect when a file is malformed?"""
file_name = 'resources/sample_malformed.fits'
file_path = pkg_resources.resource_filename('pypahdb', file_name)
file_path = importlib_resources.files('pypahdb') / file_name

self.assertRaises(OSError, Observation, file_path)

Expand Down

0 comments on commit ebcdbba

Please sign in to comment.