minor fixes

mocksipipeline/util.py

@@ -8,7 +8,7 @@
 import numpy as np
 import xarray
 from astropy.wcs import WCS
-from ndcube import NDCube
+from ndcube import NDCollection, NDCube
 from ndcube.extra_coords import QuantityTableCoordinate
 from overlappy.io import read_overlappogram
 from overlappy.util import strided_array
@@ -20,6 +20,7 @@
     'read_cube_with_xarray',
     'write_cube_with_xarray',
     'dem_table_to_ndcube',
+    'build_moxsi_collection',
 ]
 
 
@@ -166,3 +167,26 @@ def dem_table_to_ndcube(dem_table):
                                         names='temperature',
                                         physical_types='phys.temperature')
     return NDCube(em, wcs=tab_coord.wcs, meta=dem_table.meta)
+
+
+def build_moxsi_collection(results_dir, all_components_sum=True):
+    """
+    Build a collection of cubes holding different MOXSI components from a directory of FITS files
+
+    Parameters
+    ----------
+    results_dir: path-like
+    all_components_sum: `bool`
+        If True, it is assumed there is a file called "all_components.fits" in this directory
+        and every member of the collection will use that data array, just with a different WCS.
+    """
+    results_dir = pathlib.Path(results_dir)
+    image_dict = {f.stem: read_overlappogram(f) for f in results_dir.glob('*.fits')}
+    if all_components_sum:
+        all_components = image_dict.pop('all_components')
+        return NDCollection(
+            {k: NDCube(all_components.data, wcs=v.wcs, meta=v.meta, unit=v.unit)
+             for k,v in image_dict.items()}
+        )
+    else:
+        return NDCollection(image_dict)

mocksipipeline/visualization.py

@@ -101,6 +101,7 @@ def plot_labeled_spectrum(spectra_components, line_list=None, source_location=No
     x_lim = kwargs.pop('x_lim', None)
     y_lim = kwargs.pop('y_lim', None)
     log_y = kwargs.pop('log_y', True)
+    ls = kwargs.pop('ls', '-')
     threshold = kwargs.pop('threshold', None)
     figsize = kwargs.pop('figsize', (20,10))
     skip_component_labels = kwargs.pop('skip_component_labels', [])
@@ -114,9 +115,9 @@ def plot_labeled_spectrum(spectra_components, line_list=None, source_location=No
         label = labels[k] if labels else k
         color = colors[k] if colors else None
         if ax is None:
-            ax = v[0,0,:].plot(label=label, color=color)
+            ax = v[0,0,:].plot(label=label, color=color, ls=ls)
         else:
-            v[0,0,:].plot(axes=ax, label=label, color=color)
+            v[0,0,:].plot(axes=ax, label=label, color=color, ls=ls)
         if line_list and source_location and k not in skip_component_labels:
             annotate_lines(ax,
                            v[0,0,:],
@@ -153,14 +154,16 @@ def plot_detector_image(moxsi_collection,
     figsize = kwargs.get('figsize', (20,20))
     cmap = kwargs.get('cmap', 'hinodexrt')
     norm = kwargs.get('norm')
+    wave_index = kwargs.get('wave_index', 0)
+    draw_prime_key_limb = kwargs.get('draw_prime_key_limb', False)
     prime_component = moxsi_collection[prime_key]
     # Build figure
     fig = plt.figure(figsize=figsize, layout='tight')
     ax = fig.add_subplot(projection=prime_component[0,...].wcs)
     if norm is None:
         _, vmax = AsymmetricPercentileInterval(1, 99.9).get_limits(prime_component.data[0,...])
         norm = ImageNormalize(vmin=0, vmax=vmax, stretch=LogStretch())
-    prime_component[0,...].plot(axes=ax, cmap=cmap, interpolation='none', norm=norm)
+    prime_component[wave_index,...].plot(axes=ax, cmap=cmap, interpolation='none', norm=norm)
 
     # Ticks and direction annotationes
     grid_color = kwargs.get('grid_color', 'w')
@@ -212,26 +215,27 @@ def plot_detector_image(moxsi_collection,
             'spectrogram_slot_0',
             'spectrogram_pinhole_0',
         ]
+        if draw_prime_key_limb and prime_key not in zero_order_components:
+            zero_order_components += [prime_key]
         limb_color = kwargs.get('limb_color', 'w')
         for k in zero_order_components:
-            _wcs = moxsi_collection[k][0,...].wcs
+            _wcs = moxsi_collection[k][wave_index,...].wcs
             px, py = _wcs.world_to_pixel(limb_coord)
             ax.plot(px, py, ls='--', color=limb_color, lw=0.5)
 
     # Add wavelength annotations
-    annotate_kw = {
-        'textcoords': 'offset points',
-        'color': 'w',
-        'arrowprops': dict(color='w', arrowstyle='-|>', lw=1),
-        'horizontalalignment':'center',
-        'verticalalignment':'center',
-        'rotation':90,
-        'fontsize': plt.rcParams['xtick.labelsize']*0.8,
-        'weight': 'bold',
-    }
-    annotate_kw.update(kwargs.get('annotate_kwargs', {}))
-
     if line_list:
+        annotate_kw = {
+            'textcoords': 'offset points',
+            'color': 'w',
+            'arrowprops': dict(color='w', arrowstyle='-|>', lw=1),
+            'horizontalalignment':'center',
+            'verticalalignment':'center',
+            'rotation':90,
+            'fontsize': plt.rcParams['xtick.labelsize']*0.8,
+            'weight': 'bold',
+        }
+        annotate_kw.update(kwargs.get('annotate_kwargs', {}))
         annot_pt = moxsi_collection['filtergram_1_0'][0,...].wcs.array_index_to_world(
             *np.unravel_index(moxsi_collection['filtergram_1_0'].data[0].argmax(),
                             moxsi_collection['filtergram_1_0'].data[0].shape)

pipeline/config/config.yml

@@ -12,3 +12,4 @@ apply_electron_conversion: True
 apply_gain_conversion: True
 include_psf: True
 include_charge_spreading: True
+max_spectral_order: 100  # High enough to include all orders by default

pipeline/workflow/Snakefile

@@ -9,7 +9,8 @@ instrument_design = getattr(mocksipipeline.instrument.configuration, config["ins
 
 XRT_FILTERS = ['Be-thin:Open',]
 AIA_FILTERS = ['94', '131', '171', '193', '211', '335']
-CHANNEL_NAMES = [f'{c.name}_{c.spectral_order}' for c in instrument_design.channel_list]
+CHANNEL_NAMES = [f'{c.name}_{c.spectral_order}' for c in instrument_design.channel_list
+                 if abs(c.spectral_order)<=int(config["max_spectral_order"])]
 
 
 rule sum_components:

pipeline/workflow/run_pipelines.sh

@@ -0,0 +1,24 @@
+#!/bin/zsh
+# A place to put multiple pipeline commands if you need to run one after another
+
+for element in fe si mg ne o
+do
+    # Coronal
+    snakemake ../results/cdr/rfa/ar_1h_coronal_${element}/detector_images/all_components.fits \
+              --config root_dir=../results/cdr/rfa \
+                       output_dir=ar_1h_coronal_${element} \
+                       instrument_design=moxsi_cdr \
+                       spectral_table=../../mocksipipeline/spectral/data/chianti-spectrum-${element}-sun_coronal_1992_feldman_ext.asdf \
+                       max_spectral_order=3 \
+              --cores 1 \
+              --rerun-triggers mtime
+    # Photospheric
+    snakemake ../results/cdr/rfa/ar_1h_photospheric_${element}/detector_images/all_components.fits \
+              --config root_dir=../results/cdr/rfa \
+                       output_dir=ar_1h_photospheric_${element} \
+                       instrument_design=moxsi_cdr \
+                       spectral_table=../../mocksipipeline/spectral/data/chianti-spectrum-${element}-sun_photospheric_2015_scott.asdf \
+                       max_spectral_order=3 \
+              --cores 1 \
+              --rerun-triggers mtime
+done

pipeline/workflow/scripts/calculate_dem.py

@@ -86,20 +86,15 @@ def calculate_response_kernels(collection, temperature, spectral_table):
             wavelength = trf.channel_wavelength.to_value('ph Angstrom') * u.angstrom
             gain = wavelength.to('eV', equivalencies=u.equivalencies.spectral()) / u.photon
             gain /= (trf.ev_per_electron * trf.ccd_gain_right)
-            # NOTE: the xrtpy package makes use of the new unit in astropy DN which is much more commonly used
-            # in solar physics. However, DN is not a unit recognized in the FITS standard so we substitute it
-            # for count. This is currently the unit that we use in sunpy in place of DN.
-            gain = gain * u.count / u.DN
             response = ea * gain * pix_solid_angle
             response *= get_cross_calibration_factor(key)
-
         else:
             raise KeyError(f'Unrecognized key {key}. Should be an AIA channel or XRT filter wheel combination.')
         T, tresp = compute_temperature_response(spectral_table, wavelength, response, return_temperature=True)
         kernels[key] = np.interp(temperature, T, tresp)
         # NOTE: This explicit unit conversion is just to ensure there are no units issues when doing the inversion
         # (since units are stripped off in the actual calculation).
-        kernels[key] = kernels[key].to('cm5 ct pix-1 s-1')
+        kernels[key] = kernels[key].to('cm5 DN pix-1 s-1')
 
     return kernels
 
@@ -110,9 +105,9 @@ def _model(self, alpha=1.0, increase_alpha=1.5, max_iterations=10, guess=None, u
         errors = np.array([self.data[k].uncertainty.array.squeeze() for k in self._keys]).T
         from demregpy import dn2dem
         dem, edem, elogt, chisq, dn_reg = dn2dem(
-            self.data_matrix.value.T,
+            self.data_matrix.T,
             errors,
-            self.kernel_matrix.value.T,
+            self.kernel_matrix.T,
             np.log10(self.kernel_temperatures.to_value(u.K)),
             self.temperature_bin_edges.to_value(u.K),
             max_iter=max_iterations,
@@ -122,7 +117,8 @@ def _model(self, alpha=1.0, increase_alpha=1.5, max_iterations=10, guess=None, u
             gloci=use_em_loci,
             **kwargs,
         )
-        dem_unit = self.data_matrix.unit / self.kernel_matrix.unit / self.temperature_bin_edges.unit
+        _key = self._keys[0]
+        dem_unit = self.data[_key].unit / self.kernel[_key].unit / self.temperature_bin_edges.unit
         uncertainty = edem.T * dem_unit
         em = (dem * np.diff(self.temperature_bin_edges)).T * dem_unit
         dem = dem.T * dem_unit
@@ -158,7 +154,7 @@ def compute_em(collection, kernels, temperature_bin_edges, kernel_temperatures):
     dem_res = dem_model.fit(**dem_settings)
     # NOTE: not clear why there are negative values when resulting DEM
     # should be strictly positive
-    dem_data = np.where(dem_res['em'].data < 0.0, 0.0, dem_res['em'].data)
+    dem_data = np.where(dem_res['em'].data<0.0, 0.0, dem_res['em'].data)
     return ndcube.NDCube(dem_data,
                          wcs=dem_res['em'].wcs,
                          meta=dem_res['em'].meta,
@@ -174,9 +170,9 @@ def compute_em(collection, kernels, temperature_bin_edges, kernel_temperatures):
     delta_log_t = float(snakemake.config['delta_log_t'])
     temperature_bin_edges = 10**np.arange(
         float(snakemake.config['log_t_left_edge']),
-        float(snakemake.config['log_t_right_edge']) + delta_log_t,
+        float(snakemake.config['log_t_right_edge'])+delta_log_t,
         delta_log_t,
-    ) * u.K
+    )*u.K
     # Read in spectral table
     spectral_table_name = snakemake.config['spectral_table']
     if pathlib.Path(spectral_table_name).is_file():
@@ -185,7 +181,7 @@ def compute_em(collection, kernels, temperature_bin_edges, kernel_temperatures):
     else:
         spectral_table = get_spectral_tables()[spectral_table_name]
     # Compute temperature response functions
-    temperature_kernel = 10**np.arange(5, 8, 0.05) * u.K
+    temperature_kernel = 10**np.arange(5, 8, 0.05)*u.K
     kernels = calculate_response_kernels(collection,
                                          temperature_kernel,
                                          spectral_table)